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Alkusanat 



Tama diplomityo on tehty paaasiassa viimeisen vuoden aikana tyoskennellessani tut- 
kimusapulaisena Teknillisen Korkeakoulun Kylmalaboratorion teoriaryhmassa. Erityisesti 
haluan kiittaa ohjaajaani Dosentti Erkki Thunebergia, jonka tuella ja avustuksella tyo on 
viimein saatu paatokseen. Professor! Mikko Paalaselle seka laboratorion muulle henkilokun- 
nalle kiitos kuuluu inspiroivan tyoilmapiirin luomisesta ja valvojalleni Professori Martti 
Salomaalle tyon lukemisesta ja tarkistamisesta. 

Samassa tyohuoneessa kanssani pitkaan tyoskennelleille ryhmatovereilleni Juha Kopulle 
ja Risto Hanniselle annan kiitokset muun muassa mielenkiintoisista kahvipoytakeskuste- 
luista. Riston kokemus kvasiklassisen teorian kaytiinnon ongelmista oli myos arvokasta. 
Tohtori Adriaan Schakel avusti joissakin yleisissa teoreettisissa ongelmissa viimeisten kuu- 
kausien aikana. 

Olen oppinut projektin aikana paljon uusia asioita ja olen vakuuttunut, etta minulla on 
kayttoa oppimalleni jatkossakin. Kaytannossa kaiken alia kuvatusta numeerisesta tyosta 
tein itse, mutta esimerkiksi joissakin vaikeimmissa kvasiklassisen teorian analyyttisissa 
laskuissa sain paljon apua Erkilta. Tahdon korostaa hanen tarkeaa osuuttaan tyossa. Kiitan 
myos CSC - Tieteellinen Laskenta Oy:ta tehokkaiden tietokoneiden ja hyvien ohjelmistojen 
tarjoamisesta vaativimpien nmneeristen ongelmien ratkaisemiseksi. 

Lopuksi haluan kiittaa vanhempiani kaikesta siita tuesta ja kannustuksesta, jonka heilta 
olen saanut koulunkayntini ja opintojeni aikana. Lausun kiitoksen myos kaikille ystavilleni, 
jotka ovat pitaneet huolen siita, etta olen kayttanyt osan ajastani myos muihin harrasteisiin. 

Otaniemi, 17. huhtikuuta 2000 
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1 Josephson effect 

The effects related to a weak coupling between two phase-coherent quantum systems are 
often collectively termed Josephson phenomena. The name dates back to Josephson, who 
first predicted such effects (tunneling of Cooper pairs) to take place in a weak connection 
between two superconductors, a so-called Josephson junction [||]. In this Master's Thesis I 
am concerned with superfluids rather than superconductors, but the following main results 
are at least approximately true for all cases. The name "Josephson junction" is usually 
reserved only for the pure tunneling junctions in superconductors, and microbridges or 
superfluid junctions etc. are put under the more general cathegory of "weak links" Q. 



1.1 Basic concepts of superfluid weak links 

Imagine a container of superfluid, divided into two parts (1 and 2) by a thin membrane. The 
two condensates are described by some macroscopic wavefunctions, or order parameters, 
which have well-defined but different phases and (j)2- We assume that there are no other 
relevant degrees of freedom. The two sides are then weakly coupled by introducing some 
(sufficiently) small orifice(s) in the dividing membrane. If we denote the phase difference 
between the condensate wavefunctions hy (p = (f)i — (f)2, then there will be a supercurrent 
flowing through the weak link thus formed, given by the simple formula 

/(</.) =/e sin,/., (1) 

where Ic is a critical current specific to the junction. If we apply a finite pressure difference 
AP between the two sides keeping the temperature constant, there will also be a difference 
in the chemical potentials A/i = f AP, v being the specific volume V/N . The well-known 
Josephson- Anderson phase-evolution equation ||3| 

dt~ h 

then tells that (j) will change in time. As a result of this and the periodic form of Eq. ([l|), 
a constant A// results in an oscillating supercurrent through the weak link. 

The current-phase relation (|l|) is accurate only for pure tunneling junctions. For micro- 
bridge-type weak links, to which we include the superfluid ones, there are deviations from 
the sine form |Q, The general trend is that at lower temperatures, or stronger cou- 
pling, the sine will become slanted. The length scale determining the effective size and 
thus the strength of the coupling is always the supefluid coherence length, which grows 
with temperature and diverges at the transition. For large enough effective aperture sizes, 
the current-phase relation will become hysteretic (multivalued). A hysteretic /((/>) implies 
that there is dissipation taking place at a constant rate for a constant A^, since there 
then exists a net dc current component. The dissipated energy per period of 1(0) equals 
h Jq^ max{I{(j))}d(j), taking !{(()) to be the particle current. This holds because the net work 
done by the pressure difference AP in a time dt in pushing the fluid through the orifice is 

dP = .APdAT = -MN^ = -n^^dt = him<p, (3) 

dt dt dt ^' ^ ' 
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where we used Eq. (|2|) and defined I{<j){t)) = —dN{t)/dt, see Ref. ||5|. This thermodynamic 
argument thus defines the free energy of the weak link F{(j)) and relates it to the particle 
current through the equation 

The dissipation is generally explained with the concept of phase slips, introduded by An- 
derson m. A phase slip occurs when a vortex nucleates at one wall inside the orifice, crosses 
the flow under the influence of the Magnus force, and then vanishes on the other side. As a 
result, the phase difference changes (or "slips") by 27r and some energy is withdrawn from 
the flow. This happens because it costs energy to nucleate the vortex, and when it vanishes, 
its energy is dissipated as heat. The chemical potential difference is then maintained by a 
constant rate of vortex motion across the orifice: we have 



where {dn/dt)t is the average rate of phase slip events. In what follows, we shall not be 
dealing with dynamical effects of this kind. In the later sections we shall also specialise to 
weak links in superfluid ^He and concentrate on some novel peculiarities in the current- 
phase relations I{4>) found in them. 

1.2 Josephson efFect in ^He and "^He 

As soon as the superfluid phases of ^He were discovered in 1972, it seemed obvious that 
^He should display an equivalent to the Josephson-type effects found in superconductors. 
For ^He, they had already been speculated about for long |Q, but for both systems they 
kept on defying unambiguous experimental confirmation for quite some time. Evidence for 
single phase-slips taking place in ^He, which satisfied Eq. (|5|) very well, was clear by mid- 
80s ||6|. On the other hand, ideal nondissipative Josephson behavior in ^He still remains to 
be observed. 

With '^He, the earliest experiments were made at Cornell and here in Otaniemi in the 
late 70s 0. It was not until 1988 that the first succesful observation was reported from 
University of Paris in Orsay, France |^. At vapor pressure and temperatures close to Tc 
their ^He-B weak link showed nearly ideal hydrodynamic behavior, so that Eq. (jl]) seemed 
to be well satisfied. Dissipation only started to take place at lower temperatures. The basic 
reason why ideal Josephson behavior is so much more difficult to achieve in ^He than in ^He 
is in the two orders-of-magnitude difference in their coherence lengths: in '^He is on the 
order of 1.5 A whereas in ^He it is some 700 A. For ideal behavior, the dimensions of the 
weak link should be of the same size or smaller than the coherence length, and obviously 
1.5 A (atomic size) is technically very difficult to achieve. 

While many people thought the case for ^He was then settled, a group at Berkeley 
persistently kept on developing their own experiment. After another decade, in 1998, they 



reported the discovery of a new feature in the current-phase relation |12|. (The earlier 
progress of their work is described in Refs. IH, |l^, IH].) The weak link they used consisted 
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of a regular parallel array of 65 x 65 small apertures, each some 100 nm in diameter, spaced 
3 nm apart in a 50 nm thick membrane. They found that at temperatures below about 
O.STc, the weak link could be stabilised in a local minimum of energy at the phase difference 
vr. This state was seen in the I((/>)-relation as a new branch around cj) = tt, known as the 
vr branch, or the vr state. Subsequent refinements of their apparatus brought about better 
resolution, and the tt state could then be seen as a continuously evolving kink in with 



decreasing temperature ||l^. In addition, the weak link could be found in two different 
states with their own distinct current-phase relations and critical currents. Finally, it has 
recently been reported that behavior reminiscent of the tt state can be seen in a single 
narrow slit as well as the aperture array |14|. 

In the wake of this excitement, we were among the first to start doing calculations 
on the TT state. A number of previous calculations of the Josephson effect in ^He already 



existed, but none were general enough to be able to explain the observations [l^, [l^, 



m, H, H, m, H. The new attempts include Refs. |2|, |1J5|], while, a priori, the most 



convincing would seem to be Yip's work |^ and our own ||27||. This diploma thesis reviews 
our previously published results and is a continuation of that work. Among other things, 
we show that Yip's simplified model is not likely to be a correct interpretation of the vr 
state. Ours is perhaps better, but not at all flawless either. It will not always be explicitly 
mentioned, but in all our estimates and comparisons to "experiment" we will be consistenly 
assuming the parameters of the experimental aperture array of Ref. ||T^ ]. 

Much work still remains to be done, both in experiment and in theory. For example, 
most work thus far has been done on B-phase and the effects of an external magnetic field 
have not been studied in detail. The next step in increasing computational difficulty would 
be the self-consistent calculation for apertures of finite size. Furthermore, introduction of 
A-phase, or even A-B interfaces at the weak link — with or without magnetic fields — 
and so on, provide combinations for further research for years to come. For a short recent 
review on the vr state, see the article by R. Bowley in Ref. [p^] . 

1.3 Applications of superfluid weak links? 

To try and give some sort of a motivation for these investigations, I would like to note 
the following. The Josephson effect in superconductors is best known for its important 
applications in measurements of magnetic fields. Perhaps a bit surprisingly, in addition to 
its purely academic interest, the superfluid equivalent may also be of practical use in the 
future. Whereas in superconducting SQUID loops the phase is coupled to the magnetic 
field, in superfluid loops with a weak link it is, somewhat analogously, coupled to rotation. 
Based on this it is possible to use superfluid weak links for very accurate measurements of 



absolute rotation |29, pn|. The resolution obtained thus far is not quite as good as what can 



be achieved using laser interferometry (Sagnac effect), but the techniques are improving; 



see Refs. |31 
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2 Order parameter in ^He 



Let us first consider the microscopic description of the superfluid phases of ^He. This is 
important, because the "unorthodox" results to be discussed below are direct consequences 
of the more complicated order-parameter structure of ^He, as compared with, say, con- 
ventional superconductors. The discussion given here is somewhat simplified; for a more 
complete introduction, see Refs. |33, 34, 35, 36|, for example. 

Liquid ^He is a fermion system, just as electrons in a metal, with the atoms possess- 
ing a (nuclear) spin ^. Unlike electrons, the ^He atoms are neutral, but they are very 
strongly interacting. In fact, the hard-core repulsive interaction makes pair formation in 
the s-wave state impossible, so that if pairing to a superfluid state is to take place, higher 
angular-momentum states are a necessity. Whereas in metals the weak attractive interac- 
tion needed for pairing is provided by phonon exchange, in ^He it is the Van der Waals 
interaction. Strictly speaking it is not the bare ^He atoms which form the "Cooper pairs" 
in the superfluid state: the Landau "Fermi liquid" theory can be used to reformulate the 
description in terms of weakly interacting fermionic excitations, called quasiparticles \M\. 



2.1 Spin triplet and p wave 

To a good approximation, the pairing state in superfluid ^He is spin triplet and p-wave. This 
means that the Cooper pairs have the total spin s = 1 and an orbital angular momentum 
^ = 1, which is consistent with the requirement of the antisymmetry of the total fermionic 
wavefunction. Possible mixing in of other kinds of I, s combinations is usually neglected. 



2.1.1 Orbital angular momentum, I = 1 

Let {x',z = 1,2,3} denote an orthonormal triad of vectors in real space. We use these 
as the quantization axes for the orbital angular momentum: L = x'^Li -|- ^2^2 + X3L3. 
The operator I/3, for example, has the three eigenstates ji^+i), l^ii) and \Yq), with the 
eigenvalues 1, —1 and 0, respectively — the functions Y^{9, (f)) = {9, (j)\Yj:^) are the (suitably 
normalised) spherical harmonics. Prom these, we construct the following orthonormal basis 
states \Li) = -\Yl^) + \Y\), IL2) = i(|y^i) + \Y\)), IL3) = V2\Yo^), which have the 
property Li\Li) = 0,i = 1,2,3. A general normalised p-wave state may now be expressed 
as |k) = A:i|Li) -|- /c2|-^2) + ^sl-^^s), where k = ki^[ + /C2X2 + ^13X3 is a unit vector. The 
vector transformation property of the /cj's is seen as follows. 

Let the quantization axes transform under the rotation R''{n,6) as x'- =R^ ■ x^- = 
Ylii^\^\j- The basis states for this new set of axes are then obtained from the old ones 
via \L'-) = U''\Lj) = Y^- \Li)U-j, where lJ\-a,6) = exp(— iSfi • L), i.e., L is the generator 
of rotations in the angular-momentum space [|^]. Por I = 1 and this particular basis, it 
turns out that Uij = {Li\JJ^\Lj) = Rij. Rotating the vector k and the I state |k) yield 
R^--k = Zj k ■ 4 = Ei,jRijk^j and C/'|k) = kU^\Lj) = E^,J RijklU) . Thus, 
under the /-rotations C/', the components of |k) transform just like the components of the 
vector k. We may express this as C/'|k) = | i?'- k) and thus identify rotations of an / state 
with the real-space rotations of k. 
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2.1.2 Spin angular momentum, ,5 = 1 

For spin, the procedure is completely analogous. We choose another set of quantization 
axes {x^,/U = 1,2,3} and define the pair spin operator S = xfSi + + x^Sa, which 
decomposes into two parts, one for each spin: S = ^(cti + 0-2). For 5*3, we have the triplet 
eigenstates ITT), I ii) and | ti) + lit), with the eigenvalues 1, —1 and 0, respectively. 
Using these, we define \Si) ^ -|TT) + l^a) = i(|TT) + l^a) = iTi) + lit), for 

which ^^l^^) = 0, // = 1, 2, 3. A general normalised spin triplet state may be expressed 
as |d) = dil^i) + d2\S2) + d^lS^), where d = (iixf + (i2x| + d^x^ is again a vector.^] Just 
as above, rotations of the axes R^{n,9) and spin rotations U'^{n,9) = exp(— i^n • S) are 
related through IJ^\d) = \ R'^- d) and can be identified. 

2.1.3 Pair state and difTerent forms of the order parameter 

Let r be the center-of-mass coordinate of a pair. A general I = s = 1 pair state |^'(r)) may 
now be expressed in the direct-product space spanned by the above spin and orbital basis 
states: 

\P{r)) = Y,\^i^)\S,)\Li). (6) 

Here the amplitude Af^i{r), or pair wavefunction, is usually called the order parameter of 
^He. Note that we always write spin indices with Greek letters and orbital indices with 
Latin letters. Often, Eq. (^) is expressed in a spin-state form by defining |d(r,k)) = 
(k|P(r)) = k)!*?^), where = X^j^M*^* ^^'^ ^« ~ i^l^i)- If we further introduce 

Acj^(r, k) = (a/3|d(r, k)), where a, /3 =], |, and write A = [A^^], we get the most common 
form of presenting the order parameter |3^ {a^ are the Pauli matrices): 

A(.,t,=j:.„(.t)fe„.^,=(-*+* % ). (7, 

This "anomalous pair ampilitude", "off-diagonal mean field" or "gap matrix" should be 
defined more rigorously from a second-quantized viewpoint; it is a thermodynamic expec- 
tation value, and need not be interpreted as a wavefunction at all. The quantity A can 
also be written for a more general BCS type superfluid [^]. It is a second-rank spinor, and 
transforms under spin rotations as A' = UAUj' , where U_ is the 2x2 Wigner rotation ma- 
trix. Finally, I note that we can also give a "dyadic" representation for the order parameter: 
A= X^^j A^jX^x-, such that d =yl-k. This rests completely on the possibility of identifying 
the spin and orbital bases with the real space vectors defining their quantization axes, and 
therefore works only because we have I = s = 1. Since A^i transforms as a vector with 
respect to both indices, it is sometimes called a "bivector". Under simultaneous rotations 
of both the spin and orbital spaces, it transforms as a second-rank tensor. We denote the 
3x3 matrix representation of this tensor with A= [^^j] • 

^Here d is not necessarily normalised to unity, but may be proportional to an overall (complex) 
temperature-dependent "energy gap" A. Often, as below in the context of quasiclassical theory, d is written 
as A. Its name varies in the literature: gap vector, order parameter vector, spin vector, etc. 
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2.2 Broken symmetries 

We chose above the spin and orbital bases such that the corresponding amphtudes would 
transform conveniently under rotations. It can be stated that they transform according 
to three-dimensional representations of the rotation group, namely the rotation matrices. 
Apart from the rotations, another type of a (global) symmetry transformation is also 
attributed to the order parameter. This is the global gauge transformation, i.e., an overall 
shift in the complex phase of Eq. @. The two rotations plus global gauge transfomations 
form the group G = 50(3)* x 5*0(3)' x C/(l), and, neglecting any spin-orbit coupling, all 
operations of G leave the free energy invariant. But operations of 50(3)*, 50(3)' or U{1) 
separately do change ^^j, that is, the physical state: these symmetries are therefore said 
to be broken in the superfluid phases of '^He. Combined operations of G may or may not 
be symmetries of A^i, depending on the situation. For example, in the bulk B phase (see 
below), simultaneous I and s rotations remain a symmetry and form the group 50(3)'+*. 
A further example of a broken symmetry is given in the next section. 

2.3 Order parameter in the B phase 

The spin and orbital states |d) and |k) have the following properties d • S|d) = and 
k • L|k) = 0, as one can easily check. The vectors d and k thus point in directions of zero 
spin and angular momentum projections, respectively, and these directions are related by 
d =j4-k. The simplest possible case one can now consider is that where the quantization 
axes are equal, x' = x|,i = 1,2,3, and A= A / , or d = Ak, where A is (for now) a real 
valued constant of proportionality. In this case Afj_i = A(5^j and one finds that if J = L + S 
is the total angular-momentum operator, then J^|i-*(r)) = 0, i.e., the state with d = Ak 
is an eigenstate of J with the eigenvalue j = 0. 

In fact, this corresponds to an important stable bulk phase of superfluid ^He: the 
Balian-Werthamer (BW) state, or the B phase. The amplitude A = A(r) is a temperature- 
dependent "energy gap", which is related to the energy needed for breaking a Cooper pair. 
It is independent of k, which makes the B phase isotropic and in many ways similar to 
superfluid ^He, or s wave superconductors. This is in contrast to the other stable phase 
(the A phase), which is axially anisotropic, but we are not concerned with that here. 

The fact that j = is specific to the arbitrary choice that d is parallel to k. If we do a 
rotation d = A (n, 0) • k, then the resulting |-P(r)) is no longer a state of definite j, but 
rather a superposition of the j = 0, 1,2 states. However, the states obtained in this way 
are degenerate in energy with the unrotated one and should thus be equally probable to 
occur. Taking also into account an arbitrary phase factor, the most general form for the 
order parameter in the bulk B phase is A= A R {h,9)e^'^, or 

A^i{r) = AR^,{h,9)e't (8) 

A dipole-dipole interaction between the nuclear magnetic moments fixes the rotation angle 
to the value 9 = Oqk, 104°, but this still leaves the direction of the rotation axis n arbitrary. 
Near a wall, the bulk form, Eq. (|8|), is modified so that A should be replaced by a more 
general tensor quantity, as will be discussed below. 
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Figure 1: A hole of diameter D in a wall of thickness W. 



3 Ginzburg-Landau calculation 

As a first case we considered a single cylindrical aperture (see Fig. ||). The major task here 
is to calculate the order parameter self-consistently in and around the aperture. Doing this 
in the general quasiclassical formalism for all temperatures would be a tremendous task. It 
would be even more so, if the full circular symmetry of a cylindrical aperture could not be 
assumed, which will turn out to be the case. However, near Tc we may apply the Ginzburg- 
Landau (GL) expansion of the free energy and find the stationary order-parameter field 
by minimising it. This procedure has alredy been presented and described in more detail 



in Ref. |38|, but I include the results here for completeness] 



3.1 GL theory 

On small length scales (on the order of the coherence length, ^ql) one may neglect the 
weak dipole-dipole interactions between the nuclear moments. Then, in the absence of a 
magnetic field, the GL free-energy expansion includes two terms pf 



F = Fb + Fk = Jd\[h + fk], (9) 

which are called the bulk and gradient energies, respectively. The bulk, or condensation 
free-energy density is given by 

/, = -a TriAA^*) + /3i|Tr(IIT)|2 ^ p^[Tr{AA^*)f ^^^^ 
+ /^sTr (II ^ I * I ) + /?4Tr (II II ) + /^sTr (11 I * I . 



^The results presented in the special assignment were also slightly erroneous. 
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It includes all terms up to fourth order in A^i which are invariant under separate rota- 
tions of basis vectors in the spin and orbital spaces, and are linearly independent. Strictly 
speaking, in addition to the dipole energy, we are also neglecting the weak spin-orbit cou- 
pling of the Cooper pairs, both of which would violate this symmetry. The gradient, or 
kinetic energy density arises from a spatial bending of the order parameter field on the 
coherence-lenght scale: 

A = K^diA^id.A^j + K2diA*^jdiA^j + KsdiA;jdjA^i. (11) 

These terms are also rotationally invariant with respect to both indices. For a better 
justification of these forms for the energy densities, see for example Refs. |^, 

The supercurrent is related to the gradient energy by the de Gennes procedure for 
ascribing to the Cooper pairs a fictitious "gauge charge" and introducing a gauge field 
which couples to it |M, li^. For a mass current, these are the Cooper-pair mass 2m3 and 



an external velocity field Vn ||36|. To obtain the current density j, one must replace V in 
fk by the "gauge invariant" form V — i^^Vn and vary F with respect to Vn 

6F 

j(r) = -lim— (12) 

which means that j and v„ are, in some sense, conjugate variables. To the first order in Vn 
this yields for the i'th component of the mass-current density the expression 

j. = l^im [K^A;,^,A^J + K2A;^diA^j + KsA^^djA^,] , (13) 

which is exactly conserved at a minimum of Eq. (I). 



3.2 Implementation 

The minimum of the functional in Eq. (^) was computed by deriving the corresponding 
Euler-Lagrange equations, discretizing them in a 3D lattice, and solving the resulting set 
of nonlinear equations by iteration. The lattice is shown in Fig. [l|. Only a m'm2' symmetry 
(m denotes mirror reflection, 2 a rotation by 180°, and prime a time inversion]^) was used 
in the calculation. 

At the wall we required the order parameter to vanish, which corresponds approxi- 
mately to a strongly diffusely scattering surface. Although this is probably the correct 
limit in most experiments, here the choice was more a matter of convenience: any other 
boundary condition on the wall would have made the problem quite difficult. In the bulk, 
the boundary condition 

hm A^,(x,y,z) = Ae^'-^/^i (14) 

was imposed, where (j) is the phase difference between the two sides. This requires the bulk 
superfluid to be in the B phase, with the n vectors on the two sides sides being parallel. The 

^These kinds of group operations conform to the International notation of crystallographic point groups; 
they should be understood here as acting simultaneously on both the spin and orbital spaces. 
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Figure 2: The nonzero order-parameter components on the vr branch at (f> = tt. The large 
components in the middle are a product of the broken symmetry: on the normal branch 
all components would vanish there and the state would possess the full circular symmetry 
of the aperture. The wall of width W = 6^gl is shown shaded and the diameter of the 
aperture was D = IO^gl- 



general rotation matrix need not be considered here, since we may neglect the dipole-dipole 
energy on these length scales. 



The current was calculated by integrating Eq. (|13| ) over the aperture. It is important 
here to keep track of current conservation, and the self-consistency requirement that the 
total mass current through the hole should be exactly given by 

2ms dFi^ .... 

This relation should be possible to confirm at least within GL theory, and is really equiva- 
lent to current conservation in the stationary configuration which minimises the free energy. 



3.3 Results: a trapped vortex state? 

Figure ^ shows the nonvanishing components of the order parameter in a particular new 
state which we found to be stabilised for phase differences around (p = tt. The components 
Azy and Ayz, which bulge in the middle, break the circular symmetry of the aperture. The 
presence of the wall and the aperture has already reduced the bulk symmetry of the B 
phase down to where oo denotes continuous rotation (around z), and this is valid 

on any normal J((/>)-branch. But now we have a case of broken symmetry, where even the 
simultaneous spin and orbital rotations around z result in new states R^iy{z)A^jRj^{z) ^ 
A^i with degenerate energies. The remaining symmetry of the state is only m'm2\ which 
is just what we used in the calculation; assuming a symmetry somehow higher would have 
left the state undiscovered. 
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An interesting fact is that the structure of the order parameter in this state closely 
resembles that of a double-core vortex |41|. The state, which we proposed to be related to 
the experimentally observed vr state, would then correspond to a situation where a half- 
quantum vortex has crossed the aperture and formed a trapped double-core vortex state 
in it. In other words, a phase slip by vr has occurred. The analogy is perhaps more easily 
understood in the case of an infinitely long narrow slit (a purely 2D situation), where a 
similar state was found. It is partly a mystery why that was not found already in the 2D 
calculation reported in Ref. ||l9|. 

The lower part of Fig. |3| shows a phase diagram describing the stability of the new state 
with varying junction parameters, namely the thickness of the wall, W, and the diameter 
of the hole, D. In region (a), i.e., for the smallest apertures tested, the state was not found. 
In region (b) the "vr branch" appears, having a negative slope of the current-phase relation 
at (/) = TT (insert at upper left). In regions (c) and (d) the slope is positive and the vr-branch 
is very strong (insert at upper right). In (d) the state is only a local minimum of energy 
at (p = TT, whereas in (b) and (c) it is a global minimum (although just barely). Marked 
with a line are the dimensions (in units of the temperature- dependent coherence length 
^GL = / ^/WA{T)) of one experimental aperture, and the approximate observation 
point of the vr state is denoted by a cross. These are deep in the region where our vr branch 
is absent. Some of this discrepancy might be explained by the fact that the GL model we 
used works best only near Tc. But this does not seem like the whole truth: it appears that 
if this calculation is to have anything to do with the vr state of Ref. |12|, then the existence 
of a large array of holes in this experiment should play a special role. 

On the other hand, our vr state could well be what is seen in the single narrow slit of Ref. 
||l4| . However, there are some differences which require an explanation. The experimental 
J(0)'s would seem to have a vr branch even when the normal branches are apparently 
nonhysteretic. For our circular apertures and the purely 2D slit, the normal branches 
are strongly hysteretic and phase slips would appear to occur only between them. The vr 
branch in holes of this kind may thus not be experimentally attainable at all, unless one 
can somehow start directly from it at (p = vr. Smaller holes would get rid of the hystesesis, 
but in our case that would also suppress the vr state. 

Fortunately, there are also some other untested possibilities. (1) For a single slit of finite 
length, the calculation should probably be done in three dimensions — not just in the 2D 
approximation, which assumes an infinite slit length (or a channel which restricts the flow 
to be strictly two-dimensional). This is suggested by the fact that for such a slit the current 
density should be much larger at the ends of the slit than in the middle. The results could 
differ from the 2D approximation and the 3D circular hole, which are the two cases we have 
investigated so far. (2) A more general quasiclassical calculation at low temperatures could 
give a vr branch even for circular apertures with nonhysteretic normal-branch behavior. (3) 
In the experimental cell of Ref. |14|, there is no reason why the spin-rotation axes fi on 
the different sides of the aperture should be parallel to each other, which is what we have 
assumed. Instead, they may be free to point in just about any directions. 

It would be interesting to see if some of these improvements would allow direct tran- 
sitions to a vr branch (like in those reported in Ref. [19| for antiparallel n's), or even a 
continuously evolving kink. 
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Figure 3: Theoretical phase diagram for the vr state in a single aperture. The current and 
the energy as a function of phase difference (j) are also shown in two cases. Here D is the 
diameter of the aperture and W the thickness of the wall. The vr branch is found in regions 
(b)-(d), where it is locally stable at a fixed phase difference of (/> vr. In regions (c) and 
(d), it is a local minimum of energy with respect to (/> at = vr. In regions (b) and (c), it 
is also the absolute minimum of energy at (/) = vr. The parameters for one aperture of the 
experimental aperture array are indicated by the dashed line, and the observation of the 
vr state is marked by the cross. 
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4 Tunneling model 

Next we deal with the details of the "tunneling model" of Ref. j2^]. There we modeled the 
aperture array by the following phenomenological expression for its coupling energy 

Fj = -Re J;M^,M^, + + A^^;A^y)], (16) 

with a, b constants having the same phase and assumed to be real. This is obtained from 
the most general form of such coupling 

Fj = ^ a^uijAj;*A^j + c.c, {a^^ij real) (17) 

fiuij 

by introducing the restrictions that (a) the barrier be "achiral" (only states with the same 
angular momentum projections on the two sides are coupled), (b) no spin flips are produced 
(only states with same spin projection on the two sides are coupled) and (c) the barrier is 
isotropic in its plane p9[] . 

In this simple model, which I call the "18-state model" (a generalization of the Feynman 
two-state model ||5|, ^), the general superposition state describing the two superfluids can 
be expressed as 

\A) = Y,Aj:,\^,,i)L + Y,^i^^\^'^i)R' (18) 

tJ,,i tJ,,i 

where the amplitudes are A^i = A^i/A = {ij,,i\A). The basis states are just the direct 
products of the spin and orbital angular momentum states: = lid)'"^ = 

We naturally assume them to be orthonormal, i.e. — ^fiu^ij 

on both sides. In this 

basis, the tunneling Hamiltonian can be written 

Ht = -A^^[a(|/i,2:)/jL(z,/i|) +b{\fi,x)RL{x,fi\ + \fi,y)RL{y, + h.c. (19) 



Its expectation value in the state of Eq. (|l^), namely {A\TIt\A), is now equal to Eq. (|l6|). 
Note that the symmetry properties mentioned above are satisfied by Eq. ([T9|): it only 
couples states with equal I and s projections on the two sides. 

4.1 Interactions in the absence of a magnetic field 



For superfluid B phase on both sides of the barrier the coupling energy, Eq. (|16D , becomes 



Fj = - [aRj:,Rj^, + (3{Rj:,Rj!, + RjlyRj^y)] cos (20) 

Here, and henceforth, we follow the usual summation convention for repeated indices, 
except for x, y, and z. Formally this expression is obtained by substituting the order 
parameters A^-^ = AR^-^ ex.p {i(f>^'^) into Eq. (|l^), but actually there is no simple relation 
between the coupling constants a, b and a, j3. This is because the order parameter of a p- 
wave superfluid is strongly suppressed near walls and the exact meaning of Eq. (|l6|) is 
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not well defined. This will be considered in greater detail below, when we present the self- 
consistent calculations of the order parameter and a "pinhole" junction. The result of the 
quasiclassical calculation is that, near T^, the coupling energy of a dense, coherent array 
of such pinholes is of the form shown in Eq. ( ^0|) and the parameters a and /3 can be 
evaluated. More precisely, what comes out of the calculation is the the total mass current 
J, but that should be related to Fj through 

h del) ' 

as already mentioned. 

To analyse the coupling energy in Eq. (^) we have to parametrise the rotation matrices 
somehow. Typically one writes them in terms of the componenets of n and the rotation 
angle 9 

Rij{ii, 6) = cos 65ij + (1 — cos 9)hihj — smOeij^fik 

1 r- (22) 



The latter form follows, because the long-range dipole-dipole interaction |^| 

Fd = SgaA^ J dV Q + cos 9^ (23) 

fixes 9 to its minimising value 9 = 9q = arccos(— |) 104°. Other interactions are here 
not strong enough to deflect the angle from this value considerably. Even at the wall it 
costs less energy to deflect the direction of n, instead. We may write Fj in terms of 

RjliRjii = ]^[25(n^ • n^f + 30(n^ • n^) - 7] (24) 

and 

RLRuz = — [1 - 5(cos2 7?^ + cos^ 1]^) + 25 cos^ rj^ cos^ r/^ + 

+ sin 7/^ sin?7^(25cos rj^ cosr]^ + 15) cos(x^ — X^) (2^) 
— 5vl5sinr/ sinr/ (cos// —cost] ) sin(x —X )\i 

where and have been represented in the polar (77^'^) and azimuthal (x^'^) angles, 
the polar axis being the wall normal s = z. The energy only depends on the difference 

-X 



^, as required by symmetry. 
In the presence of the wall we also have the surface dipole interaction |44| 

FsD= / d2r[64(s-n)4-62(s-n)'], (26) 

where s is the surface normal. This always tends to orient n perpendicular to the wall. In 
our case, however, the wall is locally transparent at the junction, and the resulting coupling 
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interaction Fj outweighs Fg^ by orders of magnitude. But away from the transparency, 
this is the still dominant surface interaction (in the absence of a magnetic field) and should 
really leave only two choices for the orientation of the n vectors there. It is now believed 
that this explains the experimentally observed "bistability" of Ref . In our model, FgD 
is assumed to fix the left and right n's far away from the junction perpendicular to the 
wall, either parallel or antiparellel, but otherwise it is neglected — it only appears in the 
form of a phenomenological boundary condition as explained below. 

As a result, there is a competition between the orienting effects of the surrounding 
walls and the weak link. It is mediated by the gradual bending of the left and right side n 
fields between their orientations in the bulk and at the weak link. The bending energy of 



B phase is generally of the form 



Fn 



with a surface contribution 



\G1 



dRai dRaj . dRoij dRaj 



drj dri 



L 



j2 - p dRgj 



d' 



(27) 



(2^ 



The n(r) field could be calculated exactly by a minimisation of these along with the 
coupling term, but that would be overly complicated for our purposes. Instead, the following 
simplifications were madej^ 

For the surface-gradient term, we assume the GL region, 7 = 3, and = constant, 
which gives Age = 4i^A^ = 4Ag2- For the bulk terms we assume Aqi = 2Ag2, and 
expressing the rotation matrices in terms of n yields 



Fg 



tot 



:Ag2 



(29) 



This result was also obtained in Ref. |43|. Next we assume that n varies only in one plane: 
n = sinryx + cos??z. Then we average over coefficients in front of derivatives, for example 
sin^ r/(9r/)^ — > i(Or/)^. In addition, we take the average '^iCii{dirff' — > |(X]j aj| Vr/p). 
Doing all this we get 



25 



^Gtot 



G2 



j dV|V7y|2. 



(30) 



This can be minimised by a solution of the form ri{r 
similarly on the right. The divergence of the integral at r 

Ra = ay^N/^ 



+ c/r on the left side, and 
: is cut off at the radius 
4225 is the number of apertures and a = 3 fj,m is the lattice 



where 

constant of the 2D aperture array. Doing the integrals we find the quadratic forms 



F 



Gtot 



and 



pi? 

-f^Gtot 



'loo 



(31) 



where 7 = ^\/7r/VAG2fl- They contain only the polar angles of the n's at the junction 
and are thus convenient to handle. The temperature-dependent 7 can be estimated from 

■^For the meaning of the parameters you should consult Refs. |44| [i^ j, for example. 
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Figure 4: Schematic view of n(r) near the weak hnk in the vr state. Surrounding walls align 
n perpendicular to themselves, which gives rise to two possible relative orientations of n's 
far away from the weak link: parallel (a) and antiparallel (b). 



a detailed calculation but in practise it was "fitted" to experimental data, as were also a 
and /9, the only other remaining free parameters in our simple model. The infinity angles 
r/ocJ are always assumed to be either or vr to choose either the parallel or the antiparallel 
configuration. Figure § shows a schematic view of the expected spatial variation of the n 
field around the weak link. The dashed circles represent the cutoff radius Ra used in the 
definition of the 7 parameter. 

4.2 Analysis of the tunneling model 

The goal was to find the n^'^ configurations which minimise the energy Fj + Fctot for 
each fixed phase difference 4>. The following analysis is done in two steps, first forgetting 
the gradient energy altogether. It is immediately seen that the angle </) = ^ is important, 
since cost/) changes sign there and the minimisation of the bracketed term in Fj, namely 
Ej in Eq. ( |32| ) below, changes to maximisation, or vice versa. 

4.2.1 Coupling term only 

Near Tc, the coupling energy for the array was claimed to have the form Fj = — Ej cos (j), 
where the Josephson energy Ej is not just a constant, but 

Ej = [{a - m'^^.R% + m^R^^i)]- (32) 

The current through this barrier is J = JcSmcf) where Jc = {2m3/h)Ej. We first summarize 
the stability criteria of the Josephson normal and vr "branches" by taking into account this 
energy term alone. Both the parallel and the antiparallel cases are considered. 
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Consider first cos > 0. The parallel normal state with = = z has Jc = 
{2ms/h){a + 2/3) and is locally stable for all a and /3; the antiparallel normal state with 

= —h^ = — z has Jc = {2ms/h){a — |/3) and is locally stable only for a > j3. Note 
that at a = this Jc changes sign. At cos^ = 0, both normal branches become unstable 
towards a (discontinuous) jump to the tt branch. On the tt branch, which is now stable and 
has lower energy for all a,(3 > and cos^ < 0, we have again two cases 

_ J -{2mz/n)a for a > /? , . 

" \ -{2m3/h){2(3-a) for a < /? ' ^^^^ 



The first of these is achieved by 



n^'^ = -^(TV3x^y + z) (34) 



corresponding to Rj^^Rj}^ = +1, R^yR^y = —1, and Rj^^Rj}^ = —1. The other case 



IS 



n^'^ = -^(Tx + y^V3z) (35) 

with R^^R^^ = —1, R^yR^y = —1, and Rj^^R^^ = +1. That is more or less all there is to 
be said about this case. 

4.2.2 Gradient energy included 

The only other significant energy term arises from the bending or gradient energy, which 
we consider in the simplified form 



"Gtot 



7(^^-^^)' + 7(^^-^S)', (36) 



as explained above. When this is taken into account, certain changes to the stability consid- 
erations arise. The normal branch with n'^ = n^ = z (rj^ = rj^ = 0) remains unchanged 
for all a and /3. A linear stability analysis shows that the branch is now locally stable for 

5(a + f3) cos (/> + 27 > 0, (37) 

regardless of a and /3. Thus the critical phase difference (0c G [0, tt]) is given by 

cos0r = - (38) 

If 7 = this clearly reduces to the case considered previously, 0c = 7r/2. Otherwise 0c 
is moved closer to tt. The transition from the normal to tt branch can be shown to be 
continuous (second order). 

The antiparallel normal branch with n^ = — fi''^ = — z [rj^ = vr, ry^ = 0) also remains 
the same for all a > /3. For a < |/3, the critical current of the branch still becomes negative. 
However, the situation where a < /3 is more complicated, because then the bulk direction 
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of the n vector is the only thing trying to keep the branch stable. For the antiparallel case, 
a linear stability analysis shows that the branch is locally stable if 

(25a - 5/3) cos + I67 > and 

. 39 
15(a -^)cos(?;)+ I67 > 0, 

again regardless of the values of the parameters. Here comes an important point: there are 
now two critical phases whose positions depend on the parameters. For a > (5, the second 
condition is automatically satisfied and only the (upper) critical phase difference 



(25a - 5/3) 



cos '/'Ii = -77^7— (40) 



is relevant. The transition from the normal to the vr branch is discontinuous (first order) 
For a < (3 there is also a second (lower) critical phase 



15(a-/3)' 



cos^'a^ = ' (41) 



When ^/3 < a < /? the normal branch is stable between the two critical phases. Finally, 
when a < |/3 the critical phase (pel becomes irrelevant and the normal branch is stable for 
(f) > 4>c2 all the way to </) = vr. 

On the TT branch itself, the changes introduced by a finite 7 are significant, but not 
much can be said about them with analytical considerations. Whenever we are off the 
stable normal branch, the form of the current-phase relation can only be obtained with a 
numerical minimisation of Fj + Fctot- However, usually (not always) the n vectors seem 
to be directed such that the gradient energies on both sides of the junction are equal, i.e., 

= ±n^. There is some hysteresis in J{(i)) related to the discontinuous transition in the 
antiparallel case, but this is small and can also be studied only numerically. 

The cases where n is not perpendicular to the wall around = can be quite similar to 
the usual vr branches around (p = tt. There is essentially only a phase shift of tt separating 
the appearance of their current-phase relations. This should be taken into account in 
interpreting the experimental data: the "vr state" in the antiparallel case could as well 
be located at (/> = 0. We might even see a nonzero 1 — jn^l as some kind of an "order 
parameter" and the definition of the tt state, regardless of the values which cp might have 
there. The Berkeley experiment can, indeed, only determine the phase difference modulo 
TT. This is in contrast to the toroidal cell geometry of Ref. ||l^ which can resolve absolute 
phase differences. 



4.3 Results and experimental implications 

The essential features of the current-phase and energy-phase relations, with and without 
the TT state, are shown in Fig. The ratio of the gradient-energy parameter 7 to the 

^The first one of the conditions (|3e| ) is suppressed if one requires that = ±n^ always, and one 
obtains only the lower critical phase (j}c2. On the other hand, it can be seen that this requirement is usually 
statisfied on both the normal and the tt branch. This leads to the conclusion that the transition between 
the branches at <j>ci proceeds via a route where it is not satisfied and is therefore discontinuous. 
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Figure 5: Current-phase relationships and energies for the tunnehng model. The left and 
right panels correspond to parallel and antiparallel n-vectors far away from the junction, 
respectively. The directions near the junction are depicted schematically by arrows. The 
curves correspond to different values of the gradient-energy parameter 7: ideal vr state 
(7 = 0, solid line), no n state (7 = 00, dotted line), and an intermediate case (7 = 0.245 
aJ, dashed line). The parameters a = 0.2207 aJ and /3 = 0.0347 aJ are chosen to imitate 
the the experiment at T = 0.55Tc. 



coupling parameters a and /? determines their general form. If, due to a large 7, the 
n vectors are fixed exactly parallel or antiparallel and perpendicular to the wall at the 
junction, the resulting J((/>) is sinusoidal. This situation is similar to what is found in s- 
wave supeconductors where the extra degrees of freedom related to the spin-orbit rotation 
are not present. With small enough 7, there is some critical phase on the interval [0,7r] 
above which the sinusoidal branch gives way to a lower-energy state, where the n-vectors 
near the junction have been deflected from their bulk directions. The fact that the vr state 
only appears at low temperatures is explained in this model with the different temperature 
dependencies of the parameters: a,/? oc (1 — T/Tc)^, and 7 oc (1 — T/Tc), such that 7 
dominates close to Tc. 

The parameters a and (3 were estimated using a quasiclassical pinhole model. Their 
expressions are presented below in the context of a general quasiclassical calculation. Esti- 
mation of 7 involves evaluation of the gradient-energy parameter Ag2 whose calculation is 
discussed in Ref. |44| and in Ref. |48|, where Ag2 = — (^^/4m3)/?2^'°; t)e presented 
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Figure 6: Theoretical values for the parameters a, /3 and 7 in the tunneling model. These 
were obtained from Eqs. ( 1Q3D , ( |1Q4D and (^), respectively. The following experimental 
parameters were assumed: wall thickness = 50 nm, circular hole diameter D = 100 nm, 
number of holes in the lattice N = 4225, lattice constant a = 3 /xm, and total lattice area 
S = 3.8- 10~^ m^. The values plotted here were modified in subsequent calculations to give 
a better fit to the experiments. 



here. Figure ^ shows the temperature dependence of the bare a, f3 and 7, as they come 
out from their respective equations using the experimental array parameters. As can be 
seen, there is more than an order-of-magnitude difference between 7 and the coupling pa- 
rameters a and /3. This discrepancy can probably be decreased with a more realistic way 
of estimating 7. However, at least the temperature dependencies should be correct, and 
we dealt with the magnitude differences in a simple, if not completely satisfying, way. To 
provide a more reasonable fit to the experimental findings, the bare theoretical estimates 
were just scaled with some constants. Figure |^ shows the current-phase relationships for 
the best-fit scaling constants, which are given in the caption. 

In the range of parameters which appeared to reproduce the experimental results best, 
more complicated behavior than that depicted in Figs. ^ and ^ did not occur. Unfortunately, 
the more unstable behavior described in the preceding subsection is a real possibility within 
the model. It is realised for the bare pinhole estimates, for which we usually have a < (3. 
Here we avoided this problem by multiplying a with a larger constant than /3, which is an 
arbitrary choice and is not easily justified. In the more general pinhole calculation to be 
presented below, we do not have a and (3 which could be separately adjusted; there we 
have to make do with the theory as it isj^ 

^We will see that restricting the angles of the transmitted quasiparticles will favor a over f). This 
corresponds to increasing the aspect ratio W/D of the apertures. However, this not very useful, since it 
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Figure 7: Current-phase relationships for parallel (a) and antiparallel (b) n- vectors far away 
from the junction. The different curves correspond roughly to temperatures from 0.45 Tc to 
0.85 Tc with intervals of 0.05 Tc. The parameters a and /? are calculated with the pinhole 
model and 7 is estimated as explained in the text (see also Fig. §). However, to get better 
correspondence with experiments, we have multiplied the estimated values by factors of 
7.5, 1.3 and 0.15, respectively. 



Not only the relative sizes of a and f3 pose a problem. The pinhole estimates would 
appear to give slightly too low critical currents as well. This can be explained in basically 
two ways: either the experimental holes are slightly larger than those reported, or there is 
some extra leakage present at the weak link, which the experimenters are unaware of. Since 
the time these calculations were carried out and published in Ref. |^, we have learned 
that one of these cases is indeed very possible (see Sect. 7). 

Furthermore, a slow variation of the n field at the edges of the array and a slightly 
incoherent behavior of the different apertures could easily smoothen the features of J{(p) 
around the critical phases. Here the transition points are exactly defined, and appear very 
sharp and abrupt. In addition, the finite aperture sizes probably cause some extra slanting 
of J{(j)) in the experiment. To put it more generally, the specific properties of the actual 
experimental cell affect the measurements in ways that cannot be fully taken into account 
within our simple model. All of these things together could perhaps explain why the details 
of the measured J((/>)'s cannot be exactly reproduced with it. 

In spite of these problems, the tunneling model makes at least two clear predictions 
which should yield to experimental testing. The first one concerns the linear dimension L 
of the aperture array: the coupling and gradient energies depend differently on L, namely 
a/7, /3/7 oc L. This means that the larger the array is made, the more pronounced a vr 
state should be observed. The other experimental prediction concerns magnetic fields: a 
strong enough magnetic field should lock the directions of and n^, thereby making J((/>) 

also leads to a drop in the critical currents, which is not desirable. 
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strictly sine-like, although perhaps slanted or hysteretic, as is usual for finite aperture sizes 
and low temperatures. In the next subsection, the effects of an external magnetic field are 
given a bit more complete and quantitative analysis. 



4.4 EfFect of magnetic field 

In the bulk superfluid, a magnetic field will try to orient the fi vector parallel or antiparallel 
to itself. This gives rise to an energy term of the form 



(42) 



This results from a combination of the depairing effects of the magnetic field and the dipole 
energy (Hsl. A surface, on the other hand, will tend to orient n parallel or antiparallel to 



its normal, because of the surface dipole interaction in Eq. (|26D. But in a magnetic field 
there emerges another surface energy which has the form 

FsH = -d [ d?r{Ho,Raih?. (43) 
Js 

For fields greater than about 5 mT (50 G), it is usually the stronger surface interaction 



P6| . It arises from depairing effects and can be understood as follows. A magnetic field 
tends to break pairs with = along the field direction, whereas the surface breaks pairs 
with m/ = along the surface normal. Both of these cost energy, and a minimum cost can 
be achieved by choosing the direction of n such that these are exactly the same pairs. In 
the minimum of Eq. (^3|) this is satisfied as i?ai(n) then rotates the surface normal to the 



magnetic field direction |43, 4C| 



Let us first consider a homogeneous magnetic field perpendicular to the wall. In 
that case, the surface field energy is minimized with n = itz. This energy is given by 

Fli = Fi^ih = ±z) = -d\ll^\'S, (44) 

where S denotes the surface area. In the normal branch this will not affect the current-phase 
relation, since there fi is aligned perpendicular to the wall, anyway. Going to the vr state, 
however, this no longer is the case: if a > /?, say, the Josephson coupling will tend to orient 
fi such that riz = n-s = ibl/\/5 (on both sides). In this case HzRzz{nz = 1/V5)sz = and 
hence F^{hz = l/\/5) = 0, which is higher than the value in Eq. (^^. This discrepancy 
in the energies may act to suppress the tt state for high enough perpendicular fields. In the 
case of a field H" which is parallel to the wall (no perperdicular component present), Fsh 
is minimised with (fi • s)^ = (fi • H)^ = 1/5 and again the minimum is 

FgH = =nH = ±l/V5) = -d\u\\\^S. (45) 

A field parallel to the wall can thus obviously affect the normal state (fi = ±z) as well as 
the vr state since F|jj(fi = itz) = 0, which is again higher than the minimum in Eq. (^5|). 
How exactly this takes place is considered below along with a general orientation of H. 
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4.4.1 Estimation of field sizes 



To get a feeling for the magnitude of the fields which can affect the current-phase behavior 
significantly, we have to insert some numerical values. The effect of a gradient energy or any 
other interaction but FgD will be neglected. Table || gives some of the required parameters 
in the Ginzburg-Landau regime. We also know that d = doigz^'^Cch) where do ~ 2.2 for 
a diffusely scattering and do ~ 4.4 for a specularly scattering wall at zero pressure. The 
GL coherence length ^ql is defined as ^gl = \/ K/a. It can be extrapolated to lower 
temperatures with ^gl{T) = hvF/\/TOA(T), if the behavior of the temperature-dependent 
gap A(r) is known. We use the values A(0.45Tc) ^ VSksTc and A(0.8rc) \/l.75kBTc. 
Above T w O.STc, the linear approximation of Table |l] is sufficient. 



We give the field estimates assuming the conditions of the experiment ||T^. The letters 
H and L refer to the high and low critical current cases in that experiment. We denote 
the Josephson energy (in absence of the kink at i;^ = vr) with Ej Fj{'k/2) — -Fj(O). 
Furthermore, we denote E^^ = Fj{tt) — Fj{0) such that -Egain ~ 2i?j — E^^ is the energy 
gained by dropping from the (imaginary) "normal branch" to the "tt branch" at (j) = n. This 
is the energy which a (possibly additional) perpendicular magnetic field would have to win 
in order to kill the vr state and make J{(p) sinusoidal. At low temperatures (T 0.45rc), 
these energies are E'/ = 2.4 • IQ-^^ J= 1.5 eV, Eg^^ = 0.8 • 10"^^ J= 0.5 eV and E^ = 
1.0 • 10-19 J= 0.6 eV, E^^-^ = 1.1 • 10-19 J= 0.7 eV. At high themperatures (T ^ O-ST^), 
the J(0) is sinusoidal in the experiments, but we may still consider it: Ej = 0.6 • IQ-^^ 
J= 0.35 eV and Ej = 0.2 • lO-^^ J= 0.12 eV. In fact, if the gradient energy and other 
interactions but the coupling are neglected, our model can allow for nonsinusoidal behavior 
at any temperature. 

Prom the theoretical pinhole calculation, one obtains an order-of-magnitude estimate 
for the coupling parameters: 

a/m2,/3/m2 oc hypkBTcNp^^ = 3.48 • IQ-^^^. (46) 

4o^ m^ 

Here the values a = 3 ^m, D = 100 nm for the hole spacing and diameter have been 
used. For an array of 65x65 holes of this kind, the area of the array is about 3.8 • 10-^ m^ 
and the order of magnitude of the coupling energies is 1.3 • 10-11 J « 8 eV. Performing 
the actual calculation we have, for T = 0A5Tc, a 0.44 • lO-i^ J and /? 0.34 • IQ-i^ 
J. Generally speaking, these values give too small coupling energies and to get the model 
working properly, we need to multiply them with factors like 7.5 and 1.3, respectively (see 
the discussion above). These give a Josephson energy Ej on the order of 2 • IQ-i^ J 1.3 



Table 1: Ginzburg-Landau parameters for the calculation of the surface magnetic interac- 
tion at zero pressure; fe^Tc « 1.38 • 10-^^ J. 



a/(i-r/r,) 

[lO^Oj-ij^-S] 


K 

[1034j-im-i] 


gz 

[10^°T-2J-im-3] 


AV(1 - T/Te) 


1.67 


41.9 


2.33 


8.75 
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eV. The order of magnitude is then about the same as in the experiment. Based on the 
GL values in Table |l| we cannot really make accurate estimates at low temperatures. The 
temperature dependence of a and /3 in our model is also accurate only neat where the 
coupling energy then drops as (1 — T/Tc)^. On the other hand, the experimental vr state 
is only visible at low temperature, so that is the temperature region we are interested in. 
Using the extrapolation of the coherence length, we assume that we can at least get the 
correct order of magnitude out. 

If we know the total energy of the junction in its different configurations, we can 
estimate the fields which give surface energies of the same order. For example, if we want 
to know the perpendicular field which will suppress the vr state by orienting n always 
perpendicular to the wall, we can do the following: For the magnetic surface energy to win 
a coupling of the order \E\^ we should have 2|F^| > \E\. (The factor of two is present 
because the wall has two sides here.) Then Eq. (|4^ will yield a "critical field" 



IH^I- / 1^1/2 I l^l/2J lO^T (47) 

which goes as ~ (1 — T/Tc)^^^ near Tc. Now E « £'gain, as explained before. We know this 
formula will give only a very rough estimate, so we just assume E « 0.5 eV~ 0.8 • 10~^^ J. 
Using the area S = 3.8 • 10~^ m^, and the gap values mentioned above, this corresponds 
to about He ^ 6.5 mT at T = 0.45rc and He « 7.7 mT at T = O.ST^ for a diffuse wall. 
Varying the parameters and using the estimated energies for the H and L cases at these 
two temperatures, values of He in the range of about 3 mT to some 15 mT (30 G to 150 
G) are obtained. In comparison, the Earth's magnetic field is about 50 fiT (0.5 G). Our 
critical field values are at least much larger this, and it thus seems that the Earth's field 
alone should not play a decisive role. 

4.4.2 General orientations 

At a very high field, the coupling energy becomes negligible and it will be mostly the 
magnetic field which orients the n vectors. This case has been analysed in Ref. |Q. For a 
given angle between the magnetig field and the surface normal, there exist several different 
possibilities for their orientation. If the coupling energy is negligible, the orientations will 
be chosen at random between the possibilities. Each of these cases may lead to a different 
critical current. If the coupling energy is also significant, jumps between these different 
orientations may occur when the phase is varied. 

It has already been shown that a large enough perpendicular field will orient n's per- 
pendicular to the wall and make the current-phase relationship sinusoidal. This is not 
necessarily so when there is a parallel field component present. Assuming now only a par- 
allel field, oriented in the direction of ity, it can be shown that the components of the 
allowed n vectors are given by 
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where the upper and lower signs correspond to each other. On the other hand, as has been 
mentioned before, the orientations preferred by the Josephson coupHng when cos (/) < are 

where the upper and lower signs refer to the n vectors on the two sides of the junction. 
These orientations are, however, not unique: the energy is degenerate with respect to any 
rotation around the z axis. If a small parallel magnetic field is switched on, the degeneracy 
is lifted. But, as can be seen from Eqs. (^) and (^9|), some of the states allowed by the 
two interactions agree exactly. Thus, a purely parallel magnetic field can in fact enhance 
the "tt type state", where the n-vectors are not perpendicular to the wall. Unfortunately, 
although this coincidence is quite amusing it is not very interesting, since (for large fields) 
it just results in another sinusoidal J{(p) and we have seen plenty of them. 

This concludes our study of the phenomenological tunneling model and the effects of an 
external magnetic field in the depth we have felt it necessary to consider so far. Subsequent 
sections are devoted to a fully quasiclassical study of a single pinhole, plus an enhancement 
of the aperture array calculation on this basis. Effects of external magnetic fields will be 
neglected. 
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5 Quasiclassical theory 



The microscopic description of normal metals and superconductors has been formulated 
completely in terms of Green's functions [^. The problem with the full interacting Green's 
function is that is contains a lot of detailed microscopic information which is often not 
needed for practical calculations. It was first shown by Eilenberger and Larkin and 
Ovchinnikov |47| that Gorkov's equations for the (stationary) full Green's function can be 
transformed to "transport-like" equations for a quasiclassical Green's function, or prop- 
agator, where all the unnecessary information concerning atomic length scales has been 
integrated away at the outset. These are the so-called Eilenberger-Larkin-Ovchinnikov 
equations. 

The same approach can be used equally well for ^He which is also a degenerate fermion 
system, or a Fermi liquid A condition for the applicability of the quasiclassical 

approach is that the characteristic length scales are much larger than the Fermi 
wavelength Xp = In/kp. Similarly, energies must be much lower than the Fermi energy 
Ep = ksTp. If time dependence were to be considered, we would also require the time scales 
to be much longer than the inverse Fermi frequency h/Ep. We are, however, only concerned 
with stationary effects. The length scale in superfluid ^He is set by the coherence length 
^0 = hvp/27rkBTc, and energies by the transition temperature Tc or the gap A k^Tc. At 
low pressures, for example, we have ^ 70 nm ^ Ap 0.8 nm and Tc ^ 3 mK <C Tp ~ 1 
K, so the requirements are well satisfied. In what follows, we use a weak-coupling form of 
the quasiclassical approach, where quasiparticle-quasiparticle scattering is neclected. We 
also restrict to the vapor pressure in all the numerical estimates (see Appendix B). 



5.1 Eilenberger equation 

The Green's functions being treated here are single-particle temperature (or Matsubara) 
Green's functions written in the "Nambu matrix" representation. They are 4x4 matrices 
in a direct-product space of particle-hole and spin spaces. We denote such Nambu matrices 
with a "breve" accent. The general stationary temperature Green's function in k-space is 
of the form G(ki, k2, 6^)- Here = TrkBTc{2m + 1) are the discrete Matsubara energies, 
which are the Fourier variables of the imaginary time parameter r. The fact that (2m -|- 1) 
assumes only odd integer values is a consequence of the underlying Fermi statistics [S^ . 

The most straightforward way to derive the Eilenberger equations is to start with 
the general equation of motion for G, the Dyson equation, and do the "left-right trick". 
This begins by writing down the left and right side Dyson equations, which are, sym- 
bohcally, G-^G = 1 and GG-^ = ij] Here G"! = ^ - S, where S is the self-energy 
and Gq ^(ki, k2, Cm) = [icm — '^3'?ki]5ki,k2 is the inverse Green's function for noninteracting 
fermions, with = ^ki — Ep and = diag(l, 1, —1, —1). Then one proceeds by rewriting 
these equations in the variables k = (ki -|- k2)/2 and q = ki — k2, assuming that q <^ kp. 
After this one multiplies them with fa from left and right, respectively, and subtracts them 

'^The products are really so-called folding products, which include integration over the "internal vertices", 
and a corresponding matrix summation over the Nambu indices. Also, in the r-representation, is a 
differential operator, acting either to the right or to the left |Q. 
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to make terms containing cancel. Finally, one takes the Fermi-surface limit for the self- 
energy, defining a(k,q,em) = a'^ik^k, q, em)T3, and integrates" the whole equation, so 
that propagators will only appear in the quasiclassical form 

5(k,q,em) = — / d4kG(k, q, e^)- (50) 

The "quasiparticle renormalisation factor" a appearing here may be chosen arbitrarily. The 
cutoff energy E^. <C -Bp is also arbitrary, but the results should not depend on how it is 
chosen: G is nonzero only near the Fermi surface. A further Fourier transformation from q 
to r plus a small rearrangement of terms puts the Eilenberger equation into the convenient 
final form 

[i<^mh - 5] + i/iVFk • Vr5 = 0. (51) 

This is a first-order differential equation for the Matsubara propagator ^(k, r, e^) along 
classical trajectories of quasiparticles. Some information concerning the normalisation of 
g is lost in the left-right subtraction process and, therefore, a separate normalisation con- 
dition has to be introduced |H- With a suitable choice of a this condition, which is to be 



satisfied by all physical solutions of (|5l[) , can be written 

gg = -i. (52) 



To give a closed system of equations ( pi] ) and (52) still need to be supplemented by some 



self-consistency equations for the self-energy a. These are considered below. 



5.2 Impurity scattering 

The quasiclassical theory is perhaps most convenient for solving static problems which 
involve nonuiformities, such as walls, interfaces or impurities. But the above derivation is 
not really valid then, since for these kinds of strong scatterers the condition g <C fcp is 
not satisfied. The proper way to proceed, then, is to utilise the formalism of scattering 
t matrices. This generally leads to a mess, and I cannot go very deep into it here. The 
general idea is that the self-energy is divided into two parts Ti + V , where S contains the 
weak interactions, like the pairing amplitude, and V is the strong scattering potential. The 
perturbation series for G can be written in the recursive form G = Gq + GoiYi + y)G, 
which is then decomposed into three parts |50| 

G = Gi + GiTGi (53) 
f = V + VGif (54) 
Gi = Go + GoSGi or Gf ^ = Gq ^ - S(G). (55) 

Here an unphysical intermediate Green's function Gi has been introduced for convenience: 
the corresponding quasiclassical gi can be solved with Eq. (|5l|), with the effect of V coming 
to play only through G in a self-consistency equation S = S(G). If (an approximation 
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for) the t matrix can be found, then g, the quasiclassical counterpart of G, can be solved 
iteratively. The t matrix Eq. (^) is essentially the Lippmann-Schwinger equation, presented 



in books of standard quantum mechanics |51|. The quasiclassical form for this equation 
will depend on the type of the scatterer; later on we are concerned with the t matrix of a 
specularly scattering wall. The quasiclassical forms for T and V are obtained by taking the 
Fermi-surface limits, t(k., k', e^) = aT'(/cpk, A^pk', )r3, and v(k,k') = aVikpk, kFk')T3, 
as in the case of the self-energy a. All propagators are transformed using Eq. (|50|). 

5.3 Self-consistency equations 

The propagator and the self-energy are usually decomposed into "scalar" and "spin vector" 
components, as can be done for an arbitrary Nambu matrix. We write 



9 + Sl2L (/ + f • ^)i£2 



a 



u + u ■ a 
ia;2 A* • a 



A + z^, 



(56) 
(57) 



where g_ = xct^^ + yoji + z£3, and CTj, i = 1, 2, 3 are the Pauli matrices. There should also be 
unit matrices \ multiplying the scalar components, but they are omitted for brevity. The 
upper right block in Eq. ( ^7|) now contains the pair amplitude from Eq. (|^), and this is 
where the information about the order parameter comes in. The self-consistency equations 
for the off-diagonal self-energy A (the gap vector A) and the diagonal self-energy v (the 



"Landau molecular fields" z^, v and v^v^ which are all real valued) are |48 

d^k' 



z.(k,r) 
i.(k,r) 
A(k,r) 



m 



47r 
47r 



A^(k-k')<7(k',r,e„ 
A«(k-k')g(k',r,e, 
y(k-k')f(k',r,e^: 



(58) 
(59) 
(60) 



with t'(k) = i^(— k) and i>'(k) = ^'(— k). A scalar order-parameter part A(k, r) is missing 
because of triplet pairing. If we were to consider scattering from impurities, an additional 
impurity self-energy p would have to be added to (|57|). In the absence of mass currents v 
vanishes, and in the absence of spin currents v vanishes. For simplicity, we assume both of 
them to equal zero, although there are always spin currents flowing along surfaces in the 
B phase, which are what we are interested in ||5^. The remaining self-consistency equation 
for the order parameter A can be put into many different useful forms, such as 



T 



Ahi— +7rA;BrJ^ 

c 



47r 



f(k',r,e™)(k'-k) 



0. 



(61) 
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5.4 General symmetries 



We denote 2x2 matrices with an underline; for example the Pauli matrices in spin space 
are ct^, i = 1, 2, 3. Similarly, the Pauli matrices in particle-hole space are denoted as r^, and 
in the 4x4 Nambu space they take the form Tj = (8) 1, i = 1, 2, 3, where denotes a 
direct product. With these notations, the propagator and the self-energy satisfy the basic 
symmetry relations (and other forms can be easily derived from these) 



[u{k, r, e^)]^* = T3u{k, r, -e„)r3 
[u(k, r, e™)]^ = T2'u(-k, r, -em)T2 



(62) 
(63) 



where u is either g or a. They follow from the operator properties defining the Nambu 
matrices G and S, see Refs. [ 



57|| . In addition, the Eq 



More precisely, if one makes the transformations 



y) possesses the symmetry 
T2u(k,r,em)T2- (64) 



±T2g'T2 and a 



T 



-r2CT'T2 in Eq. 



(51), the equation is still satisfied for g' and a' . If the solutions are to be unique, we should 
have these satisfied as symmetries: g' = g and a' = a. The minus sign corresponds to 
physical solutions and the plus sign to unphysical solutions, which also satisfy gg = 
instead of Eq. ([5^). The off-diagonal self-energy (order parameter) is not affected by this 
requirement, but it generally causes some restrictions on a. Note, for example, that any self- 
energy term proportional to 1 does not satisfy the symmetry. This is not a real restriction 
either, since such terms would cancel in the commutator of Eq. ([5l| ) anyway. If we redefine 
the propagator components by writing Eq. (|56[) as 



c + d+ (c + d) - a 
1^2(0 — 6 -f (a — b) • ct) 



(a + 6 + (a -|- b) • a)ig_2 
c — d — 0^2(0 — d) • crcj g. 



(65) 



then the basic symmetries in Eqs. (|6^) and ([63|) can be written componentwise 



a( 


-k) 


= +a(k)* 


a(-em) 


= +a{em)* 


b{ 


-k) 


= -Kk)* 


bi-em) 


= -K^m)* 


c( 


-k) 


= +c(k)* 


c(-em) 


= +c{emT 


d{ 


-k) 


= -d(k)* 


d{-em) 


= +d{emT 


a( 


-k) 


= -a(k)* 


a(-em) 


= -ha(em)* 


b( 


-k) 


= +b(k)* 


b(-e„^) 


= -b(e„)* 


c( 


-k) 


= +c(k)* 


c(-em) 


= +c(e^)* 


d( 


-k) 


= -d(k)* 


d(-e„^) 


= +d(e^)* 



(66) 

Further symmetries follow from the geometry of a specific problem through the order 
parameter and the Eilenberger equation. 



32 



5.5 Physical and unphysical solutions — the multiplication trick 



As long as the symmetry in Eq. (|6J) is satisfied by the self-energy, the redefinition in Eq. 
(65) leads to a significant simplification in the calculations. This is because Eq. (^) can 
then be decomposed into three independent blocks of equations 







duC = 





(67) 


icmb + iAi 


c + 


-hvpdua = 





(68) 


ie^a + Ar 


• c + 


^hvpOub = 





(69) 


-AR6 + iAi 


c + 


^hvpduC = 





(70) 


icmb + iAid — iApj x 


d + 


^hvpOua = 





(71) 


ie^a + Arc? + Aj x 


d + 


^hvFduh = 





(72) 


- Ar • b + iAi 


a + 


^hvpOud = 





(73) 


i Ar X a + Ai X 


b + 


^hvFdud = 


0, 


(74) 



where A = Ar + IAi and u parametrises a quasiparticle trajectory: r = ro + uk . The first 
equation forms its own block, and we may always set c = 0. The Eqs. (71)-(74) form the 
physical block whose solutions include the physical solution, which only has a, b, d and d 
nonzero. The normalisation gg = —1 for physical solutions takes the component form 



—idd + a X b = 
(i^ + dd-aa + bb = -l. 

In the bulk (constant A and g), the normalised physical solutions are 



(75) 



d= , -''^ , a= ^ , b= ^ , d = 0. (76) 

To be exact, these forms are valid only for so-called unitary states, for which AA+ oc 1, 
or A X A* = 0. Fortunately, the unitarity requirement is automatically satisfied for the 
B phase, as well as the A phase of ^He 48 1. A more compact form for Eq. ([7^ ) can be 
seen directly from Eq. ( [5l| ) by setting the gradient term to zero and solving for g: 

^ iemh - A 
\/4 + |A|2 

To satisfy gg = —1 for this, one must require AA = — |Apl and note that {A, 7-3} = 0. 

However, the unphysical block of Eqs. (|6^)-([7Q|) is even more useful in practise than 
the physical one. This is because the physical solutions of Eqs. (|7l|)-([7^ can be obtained 
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numerically more easily by finding the exploding and decaying unphysical solutions a, b 
and c and then using the so-called multiplication trick. 

If we assume the order parameter to be real, A = Ar, we may choose all of the 
propagator components to be real, except for c and d which must then be purely imaginary. 
If we also assume A to be constant, it is easy to see that the unphysical block has the 
exploding (<, upper signs) and decaying (>, lower signs) solutions 



r a(0) 1 




Cm 




5(0) 


= Cl 


T\/4 + |A|2 


exp 


c(0) 




iA 





±2^el + \A\ 



(78) 



where ci is an undetermined constant. (The unphysical block has also three constant 
solutions, but they are of no interest here.) Similar "exponential" solutions exist for a 
general A as well, although they no longer have this simple form. Such solutions always 
satisfy the normalisation condition gg = 0, or in component form —a? + 6^ + c • c = 0. The 
multiplication trick now consists of the following property: the physical solution along a 
given trajectory can be obtained by taking the commutator of the exploding and decaying 
solutions on that trajectory 



IS 



aiu) = - [g<{u),g>iu)]. 



(79) 



Here the normalisation is given by 



-^{9<iuo),g>{uo)} 
-a^a^ + 6<6> + c< • c>. 



(80) 



which is invariant along the trajectory, i.e. uq may be chosen freely. The commutator gets 
conveniently rid of all terms proportional to the unit matrix and the normalisation ensures 
that only the relative proportions of the components of the unphysical solutions make a 
difference. In component form, the commutators become 



a = is[c<6> — 6<c>] 
b = is[c<a> — a<c>] 
d = is[a<6> — 6<a>] 
d = — s c< X c>. 



One can readily check that this procedure works at least for Eqs. ( [76| ) and (78). For a more 
detailed justification of this trick see, for example, Ref. 
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6 The pinhole model 

In the pinhole model one assumes a wall separating two volumes of ^He-B, whose thickness 
W is much less than the coherence length ^o- In this wall, the junction is formed by a hole, 
whose diameter D is also ^ ^o- This situation is depicted in Fig. |8|. The experimental holes 



of Ref. |13| do not exactly count as pinholes, but are not too far from them. The reason 
for considering a pinhole here is in the relative simplicity of the resulting calculations. The 
junction can be treated as a small perturbation whose effects on the order parameter can 
be seen in the associated energies only in the order 0([open area]^) = 0{D'^). Such effects 
can be neglected to a first approximation. A pinhole was first considered as a model for 
superconducting microbridges and previous calculations for ^He also exist IS, 21].|^ 
Here we generalise the previous calulations to find the current-phase relations for both 
parallel and antiparallel n vectors at the wall. A clear presentation on the critical currents 
of the pinhole model has also been lacking for a long time, and this situation is corrected 
here. Finally, in the next section, we generalise the "tunneling model" by performing the 
corresponding calculations directly from the pinhole model for all temperatures. 



6.1 Symmetries of the problem 

The presence of a wall imposes symmetry restrictions on the order parameter in the orbital 
space. If the spin and orbital axes are first taken to be nonrotated relative to each other, 
the gap vector A(k, r) has the form A^^^^ (A_|_ and Ay are real) 

A(0)(k,z) =A^(z)zz.k + A||(z)(pp + xx)-k, (82) 

where {z,p,x} is an orthonormal triad, same for the spin and orbital spaces and with z 
perpendicular to the wall. This form is invariant under rotations of the triad around z, and 
in the bulk it should have the simple limit A^'^^(k) = Ak, where A = A_|_ = Ay = const. 

^Recently, pinhole calculations for d wave superconductors have also been published Is^, H. 
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is the temperature-dependent B-phase bulk gap. The general form of A is then obtained 
from this by taking into account general spin-orbit rotations and overall phases. We assume 
these to be different on the left (L) an right (R) sides of the thin wall, and write 

A(k,.) = |-P(^'^')^>^°^(^'^) f-^<V (83) 
[ exp(i(/.^) R •A(°)(k,z) for z > 

Note that all of these configurations are degenerate in energy only assuming that spin- 
orbit coupling is negligible, but we do assume that. The actual z dependencies of the 
order-parameter components A|j(z) and ^±{z) have to be calculated self-consistently, and 
they depend on the type of the wall. We return to this shortly. 

Assume now that a small pinhole is made to the wall, which provides coupling between 
L and R by letting quasiparticles travel from one side to another. If the origin of the 
coordinates is placed in the center of the hole, then any quasiparticle trajectory though 
the hole can be parametrised as r = uk, where k is the direction of vp. Prom the symmetry 
of the problem it follows that 

A(°)(u) = A(°)(-u), (84) 
and due to the symmetry of the unphysical block of Eqs. (|68|) -([70|) this implies 



(85) 



These can be used to find the decaying solutions on one side from the diverging ones on 
the other, which is especially useful for u = 0. Here (and henceforth) we denote by '(0)' 
the propagator solutions corresponding to A^''-', and Eqs. (85) do not hold for the general 



(0)/ 

a> [u, 




= +«?'( 








= -»?'(■ 








= +c<»'( 





order parameter. But, in fact, we only need to solve the propagator in the simpler case. 
Writing c = ilm c, Eqs. (|68|)-([70|) become (now for an arbitrary r = tq -|- nk) 

e-mb + ]^hvYdua = (86) 
ema + A(°) ■\mc + ^hvpdub = (87) 
A^^h + ^hypdulm c = 0. (88) 

From their solutions , and c^'', we get the solutions of Eqs. (]68|)-([70|) for general 



A, Eq. (p3D, by forming the following linear combinations (on either L or R) 





(0) J 
= cos q. 


) + i6|) 




■ (0) ■ 

= ia> sm (i 











sin<p 

COS0 (89) 
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From Eq. (^), we observe that c^'^) || A*^"^, and from Eqs. ( |83| ) and (89) we find that 
also c II A, since both are obtained by the same rotation. However, note carefully that 
k A, except in the bulk. If we write k = cos??z + siniJp, we instead have the relations 
A*^°) = A_L cos??z + A|| sin??p and c^^^ = c~z'^ z, + cf^ i.e., the vectors k, A^^^ and c^^^ are 
coplanar with {z,p}. Furthermore, the norm of c is always restricted by = — 6^, so 
that one of the unknowns {a^^\ b^^\ Im c'z'^ and Im c'p^) can always be found in terms of 
the others, as long as one can be sure about the sign. Eqs. ( p6| ) are thus effectively only a 
set of three linear differential equations for three real functions (a, h and Imc^, say) and the 
solution must be integrated only for all polar angles i? of k in one plane.^ For a k rotated 
from this plane around z, the components of c^'') in the fixed basis {z, p, %} are obtained by 

the same rotation: cf'\R{7,,x) ' k) = Rij{7,,x)cf^ {^)- The scalar components a and h are 
invariant under such rotations. This, as well as Eq. (p5|), is an additional symmetry which 
followed from a specific form of the order parameter^ through the Eilenberger equation; 
again, it does not hold for c in Eq. (|^), since rotations do not usually commute. Similar 
rules apply to the physical vector and scalar components. They are useful for doing some 
angular integrations over the directions k. 

This could be stated in a different way. Under rotations of k around z, the order 

parameter transforms like A^*')(k') =Rz ■ A'''''*(k), if k' =Rz ■ k. In other words, the k 
rotation around z is equivalent to a spin rotation. Other kinds of k rotations do not have 
this symmetry. But now, all spin rotations of A^''-' resulted in the same rotation of c^^^^ 
according to Eq. (|89|) . The upper left 2x2 spin block c*^*^) (k) • of the unphysical propagator 
thus tranforms as follows: 

c(o)(k).o: = cW(i?l-k')-iiV = ii'-^-Sl-cW(k')-o:' = c(o)(k')-o:', (90) 

where the final equality follows only if Rz=R ■ Thus, only a spin rotation R^ around z can 
be "undone" with a simultaneous k rotation R-=Rz- 



6.2 Propagator at the discontinuity 

At the pinhole, the order parameter "jumps" from one value to another over a vanishingly 
small distance. As long as there is only a discontinuity of this kind and no delta-function 
potentials, the physical propagator should nevertheless be continuous along a trajectory 
crossing the pinhole, since the Eilenberger equation is of first order. Let us consider trajec- 
tories r = tik, such that u = is at the center of the pinhole. We further restrict to cases 
where k is directed from left to right; for the opposite direction the roles of L and R should 
be interchanged. Since an overall phase cannot affect any physical properties of a quantum 
system, we can also safely restrict to the symmetric case (j)^ = —0/2, cp^ = -|-(/>/2 from the 
start; all results should be independent of this choice and only depend on cj)^ — (p^ = (p. 

^In practise, it is much simpler to solve for all of the unknowns separately and use the normalisation 
just to check the accuracy. Note also that if (and only if) A_l = Ay we have Cz — ccosi?, Cp — csini!?, (i.e. 
k II A II c) such that one only needs to solve a, b and c, which is trivial and has been done analytically in 
Eq. (0). 
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We define new uppercase symbols A{k) = a<^(k, n = 0), B{k) = 6<^(k, n = 0) and 
C(k) = c^^(k, n = 0). All the physical and unphysical components for a general A can 
be expressed in terms of these. Inserting the general diverging solutions ( ]89| ) to the multi- 
plication formulas (^) such that the diverging solutions (<) are taken from the left and 
decaying ones (>) from the right, and applying the extra symmetries in Eq. ( |85| ) leads to 

a(k, 0) = i s (C^ + C^) {iA sin ^(p - B cos ^(p) (91) 

b(k, 0) = i s (C^ - C^) {A cos ^0 - iB sin ^c/)) (92) 
d{k, 0) = i s [i(yl^ + B"^) sin (p - 2AB cos (p] (93) 
d(k,0) = -s X (94) 

where we also defined C ' =R -C. The normalisation constant s satisfies s(— k) = s*(k) 
and is given by 

s(k,0) = [-{A^ + B'^)cos^ + 2iABs{n^ + C^-C'^]-^. (95) 

Note that for n vectors satisfying n'^ = itfi-^ = itz these expressions still simplify consider- 
ably. To check that the k inversion symmetries in Eq. ( |6^ ) hold for Eqs. (|9l|)-(p^), one must 
notice that the sign of (/> must be reversed on interchanging L and i?, i.e., when k — k. 
These propagators can now be used to calculate physical quantities at the junction. Most 
of all, we need the mass current. 

6.3 Mass current for the pinhole 



The general quasiclassical equation for mass-current density is |48| 

— k5t(k,r, Em)- (96) 

■III 

This can be written in terms of Re d{k) alone, and for one pinhole with open area A^ the 
current J = Ao(z • j) becomes 



J = ^oSmsvpiVpTrAiBT V / d{cosd){cos'd)p{d) / —Re 

™ ^0 Jo 27r 



d{§,xA)- (97) 



Here we add an extra p{'&) to describe the possibly different transmission probabilities for 
different trajectories; unless otherwise stated, this should be set to unity in all the formulas 
where it appears. Several forms for ci(k) can be derived. With a little trigonometric trickery, 
one finds 

d(^x6) = ^-Y {B--A-)^H<P + m + 2iAB 

^ ' ^' 4 sin2 [!(</, + 66)] + B^ cos%{<p + 59)] ' ^ ' 
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where 9 = 9{'i},x) is defined by • = cos^C^. This is interesting for analytical 
considerations, but not good for numerics since it is unnecessarily complicated to calculate 
the 9 angle in general. A less sophisticated but more easily programmable form is 

Re d{i), X, cp) - - j^^2 + ^2) COS (t>-CL. C«]2 + 4A2i?2 sin2 4> • 

For parallel and perpendicular n's we have • = C2 = ^2 — i?2 and for antiparallel 
qL . qR — q2 _ ^C2, which are the two cases of special interest. In the general case, the 
following is probably the simplest form one can get: 

• = (1 + 15(fi^ • n^)){A'^ - B^) - 5[(C • h^f + (C • n^f] 

+ (25(n^ • A^) - 15)(n^ • C)(n^ • C) (100) 
+ 5^/T5[(fi^ • C) - (A^ • C)](fi^ X fi^) • C. 

Note that the numerator is just half the cp derivative of the denominator: 

Red(k,(/>) = -^^ln|s(k,<A)P. (101) 

It would seem that, apart from a constant, an analytical expression for the energy of 
the junction could be obtained by a simple integration over 0. But unless we have equal 
rotation matrices on the two sides of the junction, there should in general also be spin 



currents present, which could contribute to the energy in some unknown way |48|. For now 
we do not need the energy and will return to the problem below. 

6.3.1 Constant order parameter case 

If we assume that the order parameter is constant all the way to the wall, the current can 
be given a neat analytical expression even in the general case. Inserting Eqs. (|78|) with 
u = into Eqs. (98) and ( |97| ) and using a tabulated formula for doing the Matsubara 



summation, one finds 



1 ^ r d'^i^ / A cos( (6 + 59 A /2)\ 
J = Ao-m.vpNpTrA J2 J ^ ^^^(('^ + ^h)/^) t^nh ^ 2fcBr ) ' 

where A = |A|. Setting % = gives a formula which is essentially that derived by Kulik 
and Omel'yanchuk for superconducting microbridges The more general form is exactly 



the result obtained by Yip as an explanation for the vr state. Unfortunately, while the 
assumption of a constant order parameter is valid for superconductors, it is poor for '^He, 
as we shall see. 

6.3.2 Parameters for the tunneling model 

We are finally ready to present how the a and /3 parameters of the tunneling model were 
obtained. In the Ginzburg-Landau region, the amplitude of A should be small. But since 
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|Ap ~ |C|2 = _ ^2|^ should have - B'^\ <^ + B^. In this hmit, we can 
integrate the current in Eq. (|97|) with respect to 4> into the form of Eq. (pO|), where a and 
/? are proportional to 

a = hypMpTTkBT / d(cosi?)(cos??Xt?) V (103) 

P=lhVFNF7TkBT [ d(cOS1?)(cOS7?)j9(7?)y (104) 

Here we hace taken advantage of the transmission probability p(i?) to calculate the pa- 
rameters for a pinhole which has the same aspect ratio as the experimental apertures. If 
we assume that any trajectory hitting the wall inside the aperture gets scattered diffusely, 
i.e., into a random direction, it does not contribute to the current. Then, for a circular 
aperture of diameter D and wall thickness W, we find 

r f(7-cos7sin7) for ^ < arctanj^ ^ 
^ ^ [ for > arctan(L»/VF), ^ ^ 

where 7 = arccos(W^/L') tan This gives the probability that a quasiparticle gets trans- 
mitted given that it hit the "open area" of the hole. Thus, a and f3 here are the original 
coupling constants a and /3 per open area So of a junction: a = SoO, (3 = Sol3. Here 
So = doS, where 5 is the total area of the array and do = -kD"^ /Aa? is the fraction of open 
area per total area in one primitive lattice cell, a being the lattice constant. We used the 
values S 3.8 • 10~^ m^, 1^ = 50 nm, D = 100 nm and a = 3 /im. The propagators in 
Eqs. ( |103| ) and ( |104D were calculated only for order parameters corresponding to a fully 
diffuse wall. All the tunneling model parameters were shown in Fig. |6|. 

6.4 Boundary conditions 

The order parameter, or the z dependence of the functions Ax(z) and A||(z) near a wall, 
is needed before we can proceed to do any of the other things just described. These have to 
be calculated self-consistently, assuming some general properties for the wall and devising 
a boundary condition for the quasiclassical propagator. The "correct" approach here would 



be to think about scattering t matrices ||55|, but this would be quite complicated 



and it has been shown that much simpler models essentially reproduce the same results 



||57| , pSj p9|l . We used perhaps the simplest of all, the "randomly oriented mirror" (ROM) 
model IQ. There are severe misprints in the original publication, so it is better to express 
the whole algorithm here anew. The unit vector s is perpendicular to the wall. 

ROM algorithm for calculating the propagator 

(1) Calculate up to the wall the solutions growing exponentially toward the 
wall on all trajectories, i.e. calculate 5< for k • s < and ^> for k • s > 0. (The 
latter can actually be obtained from the former by applying the k-inversion 
symmetries in Eq. (|66|).) 
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(2) Calculate the solutions growing exponentially out of the wall with the initial 
values 

y-w^^, . ^<(^-^0-^-) , fork.s>0 

ic^o ' {^<(k',0,e^),5>(k,0,em)} 

. ^>(^^Q'^-) , fork.s<0 

k'-s>o ' i^<(k,0,em),5>(k',0,em)} 

where denotes the point at the wall. (Again, the k • s < initial condition is 
not actually needed because of symmetries.) 

(3) Evaluate the normalised physical propagator for all directions and at all 
positions from the commutator of the converging and diverging solutions, as 
explained above. 

The specular scattering limit is obtained from this by choosing 

- " - 1^' - ^' = 2s(k • s) ^^^^^ 

1 0, otherwise, 

which just requires the propagators to be continuous on mirror-reflected trajectories: 
^^(k, 0) = 5^(k, 0), where k = k — 2s(k • s). This is a simple and intuitive case, but 
not very realistic. The totally diffuse limit can be modelled by replacing 

E^k,k' — ^/ d^k'|k'.s|. (107) 

It describes a rough wall where an incoming quasiparticle can be scattered into any angle, 
irrespective of the original direction. This limit is considered to be the most realistic one in 
most cases. Fig. |9| shows examples of the resulting order parameters for these two limiting 
cases at two temperatures. The component Aj^ is more strongly suppressed than Ay, as is 
characteristic for any wall. 



6.5 Numerics 

Computation of the order parameter is based on an iterative process. Here an initial guess is 
first taken for A^'^^. Then, using the ROM boundary condition, the Eqs. ( ^6| ) are integrated 
to get the diverging (and decaying) solutions and the multiplication trick of Eq. (|8l|) is 
applied to obtain the physical solutions. Finally the self-consistency Eq. (^Tj) is used to get 
a new approximation for A^''-* and the process is repeated until convergence is obtained. 
Having a converged order parameter and the propagator components A, B and C, the 
currents, energies and other physical properties can be evaluated by simple integrations over 
quasiparticle directions. All calculations were performed in reduced units, where lengths 
appear in units of the coherence length = fiVF/2TTkBTc and energies in units of /cb^c. 
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°0 5 10 15 20 

Figure 9: Order parameter at the wall for two temperatures, T = QATc (upper bulk values) 
and T = 0.8Tc (lower bulk values). For each temperature, the lower curve is A_l(z/^o) and 
the upper Ay {z/^q). The solid lines are for a diffuse wall and the dashed lines for a specular 
one. 



In the self-consistency equation and elsewhere, numerical integrations over the polar 
angles were carried out using a 32-point (or less) Gaussian quadrature. Accordingly, the 
diverging solutions had to be calculated for up to 16 trajectories toward (k-s < 0) and away 
(k-s > 0) from the wall. For solutions diverging toward the wall, the bulk forms of Eq. ([78|) 
were used as the initial values. The initial values for solutions diverging toward the bulk 
were obtained from the ROM prescription. Only the diverging solutions for all k directions 
needed to be calculated, because after these were known, the decaying solutions could be 
found by using the symmetries in Eq. (|6^): a>^(k) = a<^(— k), &>^(k) = — 6^''(— k) and 
Im c>^(k) = — Im c^^(— k) at each point in space. All of this had to be done for a range 
of positive Matsubara energies, the negative ones being obtained through symmetries. Ten 
energies was usually enough, although a hundred would have been no problem either, 
since the calculation is not otherwise very demanding. Summing up, all that needed to be 
caculated explicitly can be described with 5<^(i?, z^c^n). For general azimuthal angles, the 
symmetries described in connection with Eq. ( [90| ) could be applied. 

The fourth-order Runge-Kutta method was used to integrate the Eilenberger equations. 
The diverging solution is easy to find because this explicit (or "forward-type") method 
should always be unstable towards finding the solution of Eq. (|8^) corresponding to the 
largest eigenvalue of its coefficient matrix. Because of this, the diverging solution is in fact 
the only solution to be found with the method. Proceeding along the trajectory the solu- 
tions will diverge approximately as (|78|). The step size hu should then be small enough to 
satisfy 2^Je^~\^^\^\^hu/hv■p <C 1 in order to have reasonable accuracy. Nevertheless, the 
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sensibility of computing such diverging solutions numerically is a bit questionable. On long 
trajectories the solutions tend to overflow on any computer, and some intermediate rescal- 
ing of the propagator has to be done. One way to circumvent the divergencies completely 
is to rescale a, b and Im c at each step, such that a = 1 everywhere: after all, only their 
relative magnitudes make a difference when the normalised physical propagator is formed 
in Eq. ([8l|). Another way would be to absorb the leading divergence into an exponential 
by defining g{u) = g'{u) exp(Gu), where G < 2^ + \^\^ jUv-p is some positive constant. 
Then one would rewrite Eqs. ( p6| ) for ^'(m), and solve them, instead. The exponential fac- 
tors cancel in Eq. (|8l| ) and the correct physical propagator should result. For some reason, 
this approach seemed to be prone to a numerical instability towards the original diver- 
gence and was not used. It is worth investigating further, perhaps with another integration 
algorithm. 



6.6 Results for a single pinhole 

Here we consider the results for calculations of the Josephson current through a single 
pinhole aperture, assuming that the surrounding walls have fixed the spin-orbit rotation 
axes n^'^ perpendicular to the wall. As discussed above, there are then two cases: = 
(parallel) and fi^ = — (antiparallel). Three different boundary conditions on the wall 
were used for both cases: (1) a case where the order parameter was assumed to be constant 
all the way to the wall, (2) a specular wall and (3) a diffuse wall. The corresponding 
current-phase relations are shown in Figs. |lO[]l^ , with the parallel case always on the left 
and the antiparallel on the right. The (mass) currents are in units of Jq = 2m^VYN-pk-BTc x 
[open area], and only phase differences in the range [0, vr] are shown due to the symmetry 
J(27r-0) = -J(0). 



6.6.1 Current-phase relations 



(1) Constant order parameter — This is the case discussed by Yip |26|, and the current- 
phase relations shown in Fig , are exactly the same as those obtained by him. They can 
be simply plotted from Eq. (|102|) provided that one knows the temperature dependence of 
the bulk gap A. These curves are for T/Tc = 0.9, 0.8, . . . , 0.1 in order of increasing critical 
current. The parallel case is well known |[l^, but the new feature found by Yip is seen in the 
antiparallel case on the right: very close to Tc the J(i;A) is sinusoidal, but at temperatures 
below about 0.5Tc a new point on [0, vr] will appear, where J = and a very strong extra 



kink in J((/)) will form around (f) = tt. This has a simple explanation in terms of Eq. ( 102 ) 
where the phase difference (j) appears only in the combination For the antiparallel h 

vectors, 9^ depends strongly on the polar angle i9 of k (cos ~ ~ T different 
quasiparticle directions contribute to the current with a different effective phase differences. 
The currents then cancel each other in such a way that a kink will appear — this strong 
calcellation is also why the critical currents are so much smaller in the antiparallel than in 
the parallel case. For a related effect in d wave superconductors, see Ref. p^] . 

(2) & (3) Self- consistent order parameter — For the more realistic surface models, the 



results are essentially different. Fig. |ll| shows J{4>)^s for a specular surface and Fig. |l^ for 
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Figure 10: Current-phase relations for a constant order parameter. 
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Figure 11: Current-phase relations for a specularly scattering wall. 
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Figure 12: Current-phase relations for a diffusely scattering wall. 
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the diffuse surface. The curves are again for T/T^ = 0.9, 0.8, ... ,0.1 in order of increasing 
critical current; for the antiparallel case also curves for T/Tc = 0.05 are shown to emphasize 
the behavior at low temperature. The parallel current-phase relations look exactly as before, 
although their critical currents are slightly reduced. But a striking difference is seen in the 
antiparallel cases: they remain sinusoidal down to very low temperatures. As can be seen 



from Eq. (|98|) , the same current-cancellation effect is still possible in principle, but here it 
is strongly reduced: an additional kink appears only at around d.2Tc. Now it also occurs 
around = 0, instead of (/> = tt. These results have been obtained for a bare pinhole without 
adjusting any effective parameters or restricting transmission angles by the probability 
p{'&). High temperatures in the diffuse case thus correspond to the tunneling model with 
a < This is why the critical current in the antiparallel case is "negative". From this 
point of view it seems that the strong rescaling of a and /3, which made it positive in the 
original tunneling model calculation, was not well justified. 

6.6.2 Critical currents 

Figure |l3| further illustrates the critical currents Jc and the possible additional extrema 
of J((/>) for a pinhole junction in all of the above cases and as a function of temperature. 
For parallel n vectors such a plot has been published in Ref. [^], but those results were 
erroneous. In this case we see that, close to Tc, Ic{T) oc (1 — T /Tc) for the constant order 
parameter and the specular surface, and for a diffuse surface Ic{T) oc (1 — T/Tcf'^ as 



expected; see Ref. ||16|. The critical current for a constant order parameter is always the 
highest and for a diffuse wall it is the lowest. For antiparallel n vectors the roles change: 
the constant order parameter case has the lowest Jc, due to the strong cancelling between 
different quasiparticle directions, but the negative extremum around (/> = vr is nearly as 
pronounced as the positive one. For the diffuse and specular surfaces, it is seen quite 
clearly that the other extrema appear only at much lower temperatures. The dotted lines 



correspond to the high-temperature approximations obtained from Eqs. ( p. 031 ) and ( 104 ) 
for a single pinhole with p(i?) = 1 in a diffuse wall, i.e. (d -|- 213) /Eq for parallel and 
(d — ^/3)/Eo for antiparallel n's. Here the energy unit is £"0 = hvpN-pk-QTc x [open area]. 
These lines follow the correct critical currents amazingly well down to temperatures less 
than T = 0.4Tc. 

In Figs. 14 and we illustrate the effect of restricting angles with a nonzero W/D in 



Eq. (|105|) . Most importantly, increasing its value results in a drastic reduction of the critical 
currents. But for the antiparallel case, there is also some subtle fine structure involved, 
and it is interesting to compare the details of the exact results and the high-temperature 
approximation. In this approximation, the tunneling model, the critical current changes 
its sign at a = |/3 where the current-phase relation J((/') = for all cj). In the exact case, 
however, this never takes place. Instead, around values of W/D where a ^ ^(5 the current- 
phase relation develops a kink, and for a short range of W/D there exist two extrema of 
J{4>) on the interval (p G [0,vr] (see Fig. 15). This is really just Yip's "vr state", which thus 
occurs at high temperatures also, but only at considerable cost in the critical currents. 



For large W/D in Fig. 14, we see that both critical currents are positive and very close 



to each other. In terms of the tunneling model, this is easy to understand. In this limit 
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Figure 13: Critical currents for different boundary conditions. Left (Parallel n's): con- 
stant order parameter (solid line), specular surface (diamonds) and diffuse surface (cir- 
cles). Right (Antiparallel n's): solid lines correspond to positive extrema of J(0) and the 
dashed lines to negative extrema; plain curves are for the constant order parameter, and 
again diamonds for a specular wall and circles for a diffuse wall. In both figures, the dotted 
lines give the estimates obtained from the "tunneling model parameters" a and 

a » /?, and the coupling energy is Fj = —aRj^^Rj}^ cos (j)- Now, for fixed = itn-^ = itz, 
we have Rj;[^Rj}^ = 1 and Fj reduces into the same simple Josephson relation for both the 
parallel and antiparallel fi vectors. 

6.6.3 Preliminary conclusions 

Based on the above results, we now see that the simple tunneling model which we initially 
considered can reproduce most of the characteristics of a pinhole junction, or a coherent 
array of such pinholes, down to very low temperatures. We may now also safely state that 
the current-cancellation mechanism of Yip is not the mechanism underlying the vr state, 
although it is in principle quite interesting. The diffusely scattering wall is likely to be the 
closest model to a real surface, and for that a purely sinusoidal behavior is found to much 
lower temperatures than O.QTc, where the Berkeley vr state is already observed. Of course, 
there is also no way Yip's model alone could explain the fact that the weak link can be 
found in two different states, both of which show an additional kink in their current-phase 
relations. Yip had to employ a magnetic field to explain this, although the fields which 
were present in the Berkeley experiment should not have been large enough to affect the 
J((/)) significantly (see the previous discussion on magnetic field strengths). 

Although we already stated that the tunneling model should be a good description of 
the pinhole array, it is still worth while to do the array calculation more accurately. Not 
least because the antiparallel vr state in the results presented above is probably wrong due 
to the arbitrary resacaling of the tunneling parameters. 
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Figure 14: Absolute values of critical currents for increasing W/D at T = OATc (see 
Fig. 15). Solid lines are exact results and the dashed lines are obtained from the high- 
temperature approximation. The small differences in details are depicted in the inset, which 
corresponds to the region where a ^ The dash-dotted lines represent the relative sizes 
of a and (3. 
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Figure 15: Effect of restricting quasiparticle transmission in the antiparallel case at T = 
0.4 Tc. The curves correspond to W/D = 0.0,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.5 from lighter 



shades to darker. The 0.0-curve (only partly visible) is the one shown in Fig. |T2|; its critical 
current is some 30 times larger than that of the 1.5-curve, and has a different sign. The 
crossing regime corresponds roughly to the tunneling parameters a « |/3, but instead of 
going through zero, J((/') exhibits the transition by forming a "vr kink". 



47 



7 Pinhole array 



The final challenge was to construct a quasiclassical free-energy functional for the pinhole 
junction. Once a suitable form was found, the goal was to repeat all the calculations 
of the tunneling model as well as possible using that. Although the new functional is 
more accurate, it is also, in some sense, more restricted, which leads to some changes in 
the results. The idea underlying the derivation is that the small perturbations caused by 
insertion or removal of a pinhole in a nontransparent wall is analogous to, say, inserting 
ionic impurities in a bulk superfluid or superconductor. These kinds of problems have 



already been succesfully solved |^, and a very similar approach was taken here. 



7.1 The pinhole free energy 

We start from the well-known expression for the energy difference between states with one 
impurity and no impurity, V being the impurity potential | |6l| , ^ 

^^tot ^ _ lT^[in(-Go 1 + S + y) - ln(-Go ' + t)]. (108) 



To eliminate the logarithm, we may apply some form of the "A-trick" [5^. We choose to 



integrate over the strength of V by making the substitution V — > AV^ and writing 

Sn'ot ^ -ilV^' ^(Gq- 1 - S - XVr'XV = \tt j'^ ^Gif,. (109) 

Here the latter equality follows from a formal application of the t matrix equation T\ = 
XV+fxGi\V imd the relation G = Gi + GifxGi. Here G = (Gq ^ -S- AF)"! gives the full 
propagator in the presence of an impurity scattering potential and Gi is an "intermediate" 
Green's function which does not include the effect of the impurity]^ The trace operation 



Tr is defined as ||49|| 

TVF(k,k',e^) = /cbT^ j ^^Tr4F(k,k, 



2k (110) 



where F is a 4 x 4 Nambu matrix and Tr4 is its trace. Equation ( |109D is now in a form 
where the propagator can be ^-integrated directly. However, to avoid a divergence in the 
Matsubara summation, we have to subtract from Eq. ( |109D the normal-state contribution 
59.^ , which is obtained by setting S = in Eq. (jl^). We define 59. = (50*°* - 5VL^ and 
transform this to the quasiclassical form 

2 " ^ 4vr io A (111) 

X Tr4[^i(k,rij„p,em)tA(k,k,em) - (k, )tT(k,k,e„)], 



^"See Sec. 5. For some further discussion on this type of operator formalism, see for example Refs. [S^, |51 
and some of the papers listed in the references. 
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where rimp is the location of the impurity and i(k, k, em) is the forward-scattering part of 



the t matrix. This formula is simpler to use than those in Refs. |50|, since it does not 
involve an integration over r. More importantly, gi is constant in the A integration. 

Now we specialise the above approach to our particular problem. The coupling energy 
we wish to know is, by definition, the difference in energies between an open pinhole and a 
blocked pinhole. It should be irrelevant how the hole is blocked as long as the transmission 
of quasiparticles is prevented. Changing the type of blockage should only change some 
constant terms in the energy, which do not depend on the phase difference or the rotation 
matrices. We might, for example, block it with a piece of specularly scattering surface, 
which corresponds to a delta- function scattering potential V6{z) in the limit V ^ oo. The 



t matrix of this type of "impurity" is of the particularly simple form |£5|, |48|| 



2vF\k-z\ XV Ao6l y 

ixik k') = ^ ^ , (112) 

2vF\k ■ z| - XV[gi{k,z = 0, e™) + ffi(k, z = 0,em)] 

where Aq is the area of the blocking piece wall with normal z (equal to the open area of one 
open pinhole), ky = k — (k • z)z denotes the parallel component of k, and k = k — 2(k- z)z. 
On inserting this into Eq. (|1 1 1]) and performing the A-integration, we find 



(113) 



F{V) = -AohvFNFirkBTy^ [ ^|k • z|Tr4 x 

m •' 

2wF|k • z| - V[5i(k, 2 = 0, em) + 51 (k, z = 0, e^,)] 



X In ■ 



2?;F|k • z| + 2iVr3Sgn(em) 



where the normal-state propagator = i-r3Sgn(em) was obtained from the bulk form, Eq. 
(77), by setting A = 0, A = 0. Finally, in the limit V ^ cxd we have 



1 ( d^k 1 

i^pinhoie = 2^°^^FiVFvrA;Br2^ / — |k-z|ln{Det4-[5i(k,0,e„) +5i(k,0,e„)]}, (114) 

where gi is the physical propagator inside the open pinhole, which has already been "solved" 
in Eqs. (|9l|)-(|9^). Here we used the general properties Trln^ = InDetA and TieiAB = 
DetA Deti? and noted that Det[iSgn(em)T'3] = 1- 



Evaluation of the determinant in Eq. (114) is straightforward in principle, but very 



complicated in practise. So, instead of actually doing that, we considered a further simpli- 
fication. Our choice to block the hole with a piece of specularly scattering wall was already 
arbitrary, so there is no reason to stick with that in case an even simpler surface is found. It 
is now possible to think of an imaginary type of surface, which retroreflects all quasiparti- 
cle directions instead of mirror-reflecting them. It would seem at least intuitively plausible 
that even this choice should lead to the correct coupling energy terms. The t-matrix for 
such a surface is obtained from Eq. ( |112| ) by replacing the mirror-reflected directions k by 
the reversed directions — k. The same replacements should then be done in the expression 



for the energy in Eq. (114) as well, but that is all. Doing this will simplify things consid- 



erably, and the determinant becomes relatively easy to calculate; the process of evaluation 
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is described in more detail in Appendix A. The end result of the calculation is that the 
coupling energy takes the simple form 

i^pinhoie = ^^Aohv^N^TTk^TY, j ^|k " z| [in | s(k, 6^) | ^ + ln(4^4)] , (115) 

where 

|s(k)r2 ^ p2^52)cos</.-C^-C^]2 + 4^2^2g-j^2^_ ^-^-^g^ 



As expected, J = jz^-o = (2?Ti.3//i)3Fpinhoie/(?</' gives exactly the mass current in Eqs. (|96| ) 
and (99) for a pinhole of open area Aq. The spin current contribution seems to have been 



absorbed into Eq. (|115| ) completely through a proper choice of the "limits of (/> integration". 
7.2 Spin current 

A mass current is associated with broken gauge symmetry C/(l) and occurs when the phase 
of the order parameter varies in space. Such a variation increases the energy and leads to 
"phase rigidity", since the phase field tends to be as uniform as possible. Similarly, due to the 
broken S0^{?>) x 5'0*(3)-symmetry, a spatial variation of the spin-orbit rotation Rfj,i{h,d) 
(in the B phase) causes spin currents [Q. This is because there is then a position-dependent 
phase difference between "spin up" and "spin down" Cooper pairs and, although the total 
mass current due to this vanishes, there can be a net transfer of angular momentum. In our 
case, there is a discontinuous jump between the rotation matrices R^l^ and spin currents 
should in general be present. 

7.2.1 Quasiclassical expression for a spin current 

In quasiclassical theory, the spin current density has the expression 

i'lpini'^) = ^VFNpTTkBT / — kfif^(k, r, e„^), (117) 



where is the 7-component of the propagator vector g |^. For a pinhole of open area 

r. 

spin 



Ao, the spin current J]r,i^ = Ao{z • Jt"") may be expressed as 



•^spin = AohvFNpTTkBTY^ I ^k • zRe d^{k,0,em) (118) 



where is the 7-component of Eq. 

Re dj = Re{-s [C^ x C^]^} = -|sp(Re s) [C^ x C^]^ 

[-{A^ + B'^)cos(t> + C^ -C^C^ X C^]^ (119) 



[(^2 + ^2) cos ^ _ • C^]2 + 4^2^2 gin2 



We may express this in terms of the rotation matrices with • C-^ = R^iR^jCiCj and 
[C^ X C-^]^ = eaf3'yRaf^R^iCkCi. There should now exist a general way of writing the spin 
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current in terms of the energy F = ^pinhole and these rotation matrices, analogously to Eq. 
(21) for the phase difference and mass current. With the help of such a relation, one would 
be able to confirm the correctness of Eq. ( 115 ), although it is otherwise not needed here. 



7.2.2 General spin current expression at a discontinuity 

Consider again the spin triplet state of a Cooper pair, now in the form |d) = {—dx+idy)\ ||) 
+{dx + idy)| ii) + dz{\ Ti) + I it))- If we apply to this state a spin rotation exp(i^2Z • S) 
around the quantization axis z, we find exp(i^j;Z • S)|d(r, k)) = {—dx + idy) ex.p{i6 z)\ IT) 
+{dx + 'idy) exp{—i9z)\ il)+dz{\ Ti) + | it))- If the rotation angle 9z possesses a nonvanishing 
gradient into some direction, then we can see that the up-spin states have a velocity 
into that direction while the down-spin states have a velocity in the oppsite direction: 
mass currents cancel, but angular momentum is being transported. This is the essence 
of the concept of spin currents. Similar considerations apply for general spin rotations 
= On. = Oxic + OyY + Ozit, and at each point in space we define a spin velocity vjpj^^ = 
{h/2m3)^6y, analogously to the B-phase superfluid velocity = (/i/2m3) V(/>, where V(/> 
is the gradient of an overall phase. Now, just as the mass current density is related to 
by j(r) = (2m3/?i)((5F/(5V(/)(r)), we also get the current density of angular momentum by 

In our case, however, we have discontinuous jumps of the phase and the spin rotation 
angles at the junction. They are of a step-function type so that [V(p]z = 4>S{z) and [V6'^]2 = 
6^S{z), where we defined the differences (p = cj)^ — (p^ and 9^ = 0^ — . The problem 
is that the latter of these bears no meaning in this global sense: we cannot add rotation 
vectors, since three-dimensional rotations do not usually commute. But for infinitesimally 
small additional relative rotations we can safely write 69^ = 59^ — 59i^ , and we expect to 
have the total mass and spin currents at the junction to be given by 



_ 2m3 dF{(l) + » ^'^^ 



h 



and J? - ^ne^e^m 



(120) 

se=o 



Assuming now that F can be written in terms of the combinations R^^R^j as above, we 
can represent such additional relative infinitesimal rotations by the expression R^^R^-Rap, 
where Rap{86) = 5a/3 + eai3'yS9-y. By using the chain rule of partial differentiation, we then 
find from Eq. (|l2|) 

77 OF d[RiiRf^{6^p + e^f,^69^)] 



(121) 
se=o 



Jlin - '"^^^^^^S^Cf^^j- ^^^^^ 



which is nothing but 



This is certainly satisfied by the expressions in Eqs. ( |115| ) and ( |118| ) and it was easy to 
guess by mere inspection before any calculation. Note that this gives again a particularly 
simple expression, if calculated for the high-temperature form, Eq. (pO|), of the coupling 
energy F. 
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7.3 Parameters and implementation 



In the case of a single pinhole we could simply assume the n-vectors to be always parallel 
or antiparallel and perpedicular to the wall, and the energy of Eq. ( [11 5] ) was not needed. 
But as we saw in the tunneling model, if we take a large array of such pinholes, the 
total coupling energy Fj can overcome the gradient energy Fctot and there will then be 
a copmpetition between these energies to find a stable equilibrium state for each given (p. 
Here we want to find the mass currents corresponding to those equilibrium states. This is 
done by minimising the total coupling energy, namely the "sum" of -^pinhole's for all holes 



in the array, plus the same model gradient energies ( pT| ) we used in the tunneling model. 
As parameters in our calculations we have the total area of the junction S and the open 
area density do- The coupling energy Fj is obtained from Fpinhole by replacing Aq with 
the total open area So = doS, which means that we assume the pinholes in the array to 
operate fully coherently. We may also study the effect of restricted angles, so the aspect 



ratio W/D in Eq. ( 105 ) is in principle third possible adjustable parameter. The value of 
the gradient-energy parameter 7 is set by the total area S and it is given in Fig. ^, but 
this can be scaled also. We are, nevertheless, in a more restricted situation now, because 
the sizes of "a and /?" cannot be adjusted freely. 

We again parametrise Fj in terms of the polar (r/) and azimuthal (x) angles of n^'^ 
with respect to some fixed axes. The polar axis z is perpendicular to the wall. Owing to 
the symmetry of the wall, the absolute azimuthal angles and x^ do not matter, only 
their difference does. We thus choose x^ = ~X and x^ = X, i-e- '^X = ~ X^, and insert 
the expressions 

= sin 7]^ cos xx — sin rj^ sin XY + cos rj^z 
n = sin r] cos xx + sin r/ sin xy + cos rj z 

into the energy, Eq. ( [115D , where C^-C^ is given by Eq. ( |100[) . The task is then to minimise 
the total energy with respect to the three real angular variables {fj^ , , x} ■ The vectors 
attain all possible values on the intervals rj^,r]^,x ^ [Oi^]i but it is easier not to restrict 
their values in the numerical algorithm. 

The numerical minimisation olF{4>) = Fj((/))-|-FGtot was carried out using the standard 
NAG Fortran library routine E04JYF. The order parameters and propagators needed for 
the calculation of Fj were obtained in the same way as in the single pinhole case. Practically 
speaking, only the minimisation routine had to be added to the program at this final stage. 



7.4 Results for the pinhole array 

Compared to the tunneling model calculation, here we chose to use slightly different pa- 
rameters for the experimental aperture array. This is because we have subsequently learned 
that, instead of cicular holes, the aperture array had actually always consisted of approxi- 
mately square holes (at least that was the form aimed at in the etching) |^. Also, instead 
of a diameter of 100 nm, these squares were roughly of size 115 nm x 115 nm. There is 
still some doubt as to whether these are the real hole dimensions, because they have been 
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0.5 1 0.5 1 



Figure 16: Current-phase and energy-phase relations for the coherent pinhole array at 
temperatures between O.STc and 0.9Tc in intervals of O.lTc. The left panels (a),(c) are 
for parallel n vectors at infinity and the right panels (b),(d) for antiparallel ones. The 
parameters used are 5 = 3.8- 10~^ m^, do = 14.69- 10~^, W/D = 0.0, and 7 was multiplied 
by a factor of 0.1. The (p sweep was done from left to right, and this is the form of the 
curves as they came out of the minimisation routine. No numerical noise was applied 
in the minimisation, and we find strongly hysteretic behavior at the transitions, which 
makes the antiparallel J{(j)) look a bit "unappealing". In the parallel case, the n vectors are 
perpendicular to the wall around = 0, but in the antiparallel case they are perpendicular 
around = vr. The two "tt states" are thus in different places. 
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Figure 17: Hysteresis associated with the jumps between the normal and vr branches at low 
temperature (T = 0.3Tc). The arrows show the routes taken when proceeding into different 
directions. How much of this is real and how much just numerical is difficult to study. All 



parameters were as in Fig. 16 



deduced after the etching to be consistent with an experiment measuring the normal-state 
flow resistance, which should depend on the open area of the holes 

Nevertheless, due to this change, the total open area has now been multiplied by a factor 
of about 1.68, which increases the critical currents quite a lot and at least solves one of the 
previous problems right away. Fig. 16 shows the results of minimisations at temperatures 
T/Tc = 0.3, 0.4, ... ,0.9 for 5 = 3.8 • 10"^ m^, do = 1-469 • lO^^ ^nd W/D = 0.0, which is 
the new "basic configuration". In addition to these values, the gradient-energy parameter 
7 obtained from Fig. ^ has been multiplied by 0.1 to get more pronounced vr states. Note, 
however, that this is now the only unjustified scaling that is being done. 

There are some clear differences between the results of Fig. 16 and those in Figs. ^ and 



13. First of all, the antiparallel "vr state", where the n vectors are not perpendicular to the 



wall, has moved from (j) = tt io (p = 27r. Moreover, the form of the corresponding J{4>) is 
now also quite different from the parallel J{4'). We know that in the experiments the low 
and high critical current J{(j))'s actually did look quite different, and it was the low critical 
current case which had the relatively more pronounced vr state. This is just what is seen 
in Fig. 16 and, therefore, at least one more problem is now corrected, although the new 
curves are not quite as "pretty" as before. I should stress that above T = 0.4rc, essentially 
the same results as here could have been obtained also by the high temperature form of 
Fj used in the tunneling model. 

As is shown more clearly in Fig. 17 for T = 0.3Tc, there is again some hysteresis 
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Figure 18: Relative angles of the fi^'^ vectors in one minimisation sweep at T = 0.5Tc for 



parallel n's at infinity. All parameters were the same as in Fig. 16. The interval with the 
strange bumps in the middle is a region where the gradient energies arising from the left and 
right sides of the junction are not equal. This is seen as an asymmetry in the magnitudes 
of the angles rj^'^. The nonzero 2x angle in the beginning follows directly from an initial 
guess, because for perpendicular n vectors the angle is actually undetermined. 



associated with the transitions. This is hard to study properly, because the analysis has to 
be done numerically. Some of the hysteresis can always be purely numerical, such that the 
minimisation routine just gets stuck in an unstable extremum of the energy. This could 
be prevented by adding some random numerical noise in the initial n configurations for 
each minimisation, instead of using the solution for the previous (j). Unfortunately, this 
might then induce some premature jumps between branches and wipe out some of the true 
hysteresis as well. I postpone any further study of this until later (in the case it should 
turn out especially interesting.) 

Figure 18 shows the minimising angles {r/^, 77^, 2%} for the curve at T = 0.5Tc for par- 
allel n vectors in Fig. 16. Previously, in connection with the tunneling model, we mentioned 
that in the vr state the n vectors are usually directed so that the gradient energies on the 
left and right sides are equal: = ±n^. Here we have a case where this is not obeyed. Up 
to about (j)/27r = 0.35, the condition \r]^\ = \r]^\ is satisfied]^, but then a smaller energy 
configuration is found, where it fails. This is seen also in the corresponding J{4>) as an 
extra bump. Here is another "inconvenient" effect which was (unintentionally) avoided in 
the tunneling model calculation, but which seems to be present for these more realistic 
choices of parameters. 

^^No need to worry about the negative polar angle: it is used here just for convenience. 
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Many other combinations for the parameters than those used in Fig. 16 were also 
tested, but none resulted in any new kind behavior: Increasing W/D decreases the critical 
currents, which is not very useful. Increasing the scaling factor of 7 over 0.2 or so begins to 
suppress the vr states too much and the discrepancy in the relative sizes of the coupling and 
gradient energies thus remains. A larger array than that in the experiments, or a better 
way of estimating 7 would still be needed to correct this. 
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8 Conclusions and discussion 



We note that our new results for the aperture array (in Sec. 7) differ somewhat from those 
published previously (Sec. 4). The main reason why the "inconveniences" related to them 
were originally avoided (by the rescaling of a and /3) was that the resulting curves looked 
too suspicious at that time. And, admittedly, we were in a hurry to publish the results. 
However, the essential physics is still the same, irrespective of whether it is right or wrong 
— and we believe it is right. I repeat the simple predictions of our model, which should be 
easy to test experimentally: a smaller array size or a strong enough magnetic field should 
suppress the vr state. If even this does not happen, then we obviously need something very 
different to explain the experiments. Such tests have not yet been carried out, but we hope 
that the situation will soon change. 

One rather trivial, and initially plausible explanation was already given in Ref. p5|| . 
There the whole vr state was explained by a simple argument where half of the holes in the 
array were assumed to be on one branch and the other half on another branch of a hysteretic 
J(0). According to this, the vr state would have resulted from the currents of different 
holes cancelling each other — much like in Yip's explanation. But this explanation was 
soon deemed unlikely, because it required the (incorrect) assumption of hysteretic "normal 
branches". Nevertheless, it is not at all clear as to why the apertures should work absolutely 
coherently. In this work we have consistently assumed it to be the case, but mainly because 
the experimenters claim so (see Ref. ||l^]). Even to themselves, it has always been a bit of 
a mystery why this holds. But perhaps some minor incoherence is just one of the many 
possible factors which make the experimental findings different from our theoretical ones. 

In summary, our current-phase relationships would now seem to mimic the experimental 
ones quite nicely, considering the crude approximations made. They fall clearly into two 
classes according to the size of the critical current, and the explanation for this is simple: 
the parallel or antiparallel orientations of the fi vectors. Both cases show an additional 
kink at low temperatures, which is stronger for the smaller Jc case. The temperature 
dependence of the form of J((/)) is also explained in a simple fashion with the different 
temperature dependencies of Fj and -Fctot- The only remaining problem is actually the 
too large magnitude of the latter (by a factor of around 10), as obtained from our estimate. 
We should try to improve the estimate, hopefully still retaining the simple quadratic form 
of Fqxqx . The problem with too small critical currents was solved by having to increase the 
size of the apertures, and now the Jc's are perhaps even too large. But our knowledge of 
the actual aperture sizes is too inaccurate to make any conclusions based on this. On the 
other hand, the current densities in a pinhole should always be larger than for any finite- 
size aperture. The next step, hopefully after correcting also the problem with the size of 
7, could thus be the self-consistent quasiclassical calculation for an aperture of finite size. 
Whether that is needed, however, still remains to be seen. 

Interest in the vr state of the Berkeley array experiment |12| was the reason why we 
started working on this subject in the first place. The pinhole array model (Sects. 4 and 
7) is now our primary candidate for an explanation of the vr state (see Fig. |l6|). But we 
should not forget the valuable "byproducts" of this project either. The vr branch found in 
the single large aperture by the GL calculation (Sec. 3) is also very interesting (see Fig. |3|). 
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And, as already mentioned, this could in fact be the proper interpretation of the results of 
Ref. 1 14 1 measured for the single narrow slit. We are eagerly waiting to hear more news on 
the progress of these experiments. 

There is also a small chance that, instead of having holes larger than they initially 
thought, the Berkeley experimentalists have some extra leakage parallel to their aperture 
array. If this were the case, also their tt state could perhaps be explained with the single 
aperture calculation. Doing a new, more sophisticated GL calculation in the future is 
therefore also a possibility — at least in the case that the aperture array model turns out 
to fail the possible tests concerning magnetic fields and array sizes. 

Furthermore, the properties of a single pinhole have now been investigated in some 
detail, taking into account different surface models (Sec. 6). Although the pinhole problem 



is already an old one, no very clear plots of the current-phase relations (Figs. I^r2|) or 
the critical currents (Fig. |l3|) have ever been published before. With the results of this 
calculation, we also consider Yip's model for the vr state disproven. 
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A Evaluation of the determinant 



The coupling energy of the open pinhole, relative to the state blocked with a specular wall 
was found to be given by 

1 I d'^k 1 

i^pinhoie = / — |k • fi| ln{Det4- [^i (k, 0, e„) + 51 (k, 0, e^)]}, (124) 

where k = k — 2(n • k)n is the reflected quasiparticle direction. The difficult part in 
evaluating this is to find the determinant of the 4x4 matrix. The expression to be calculated 
is of the form 

Det4^[5i +52], (125) 

where in our special case gi = g{k) and ^2 = 5(k). Writing the matrix sum directly in 
component form and evaluating the determinant that way is difficult, but luckily there is 
a much simpler detour. The determinant in Eq. (|125|) can be written 



1 \^ / fl \^ / 1 „ f , 1 



y (^Det4-[5i+52]J =yDet4(^-[5i+52]j = ^ ^Det4 (^-1 + -{9i,52}J , (126) 

using the general property Det^-B = DetA Det-B and the physical normalisations gigi = 
9292 = —1- Note that if gi = ^2, this equals unity: indeed, since we require a physical 
propagator to have the normalisation gg = —1, it follows that its determinant must be 
Det4^ = ±1. This is consistent with the fact that the logarithm in the junction energy 
expression should be zero if k is replaced by k (and thus ^(k) = ^(k)), i.e., if the energy 
of the transmitting junction is calculated relative to itself, and not to that of the blocked 
junction. 

In any case, the commutator {51,^2} can be simplified further without going into any 
particulal coordinate representation. One finds that it is of the form 



^{51,52} 



5 + /io:2 

/i£2 5 + i£2g • '^1^2 



(127) 



with 



g = did2 + di • d2 - ai • a2 + bi ■ b2 
g = did2 + did2 + iai X b2 - ibi x a2 
/ = di • a2 + ai • d2 + di • b2 + bi • d2 
/ = di • a2 + ai • d2 - di • b2 - bi • d2 



fl28) 



Finding the determinant from this form is easy: 



Det4(-l + ^{5i,52}) = (-//- (5-1)' + g-g) • (129) 
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In principle, it is now straightforward to just insert the expressions for /, /, g and g. In 
the present case we need the physical propagators at the pinhole, Eqs. namely 

a(k, 0) = i s (C^ + C^) {iA sin ^(1) - B cos 

b(k, 0) = i s (C^ - C^) {A cos ^0 - iB sin ^cj)) ^^^^^ 
d(k, 0) = i s [i(^2 + sin (p - 2AB cos cp] 
d(k,0) = -s X 



where 



s(k,0) = [-(^^ + 5^)cos</> + 2iABsin(/. + C^-C^]-^ (131) 



and the corresponding expressions for k. Unfortunately, the calculation turns out to be 
very tedious in practice, because the parallel and perpendicular components of C behave 
differently under the symmetry operations connecting k and k. 

However, intuition comes to rescue. As mentioned in the text, the coupling energy 
should not really depend on how the pinhole is blocked. The type of wall used to reflect 
the quasiparticles should only be seen in an additional constant in the free energy, which no 
longer depends on the phase difference or the rotation matrices. Thus, if one imagines the 
hole being blocked by some kind of material which retroreflects all quasiparticle directions, 
the interesting energy terms obtained should still be exactly the same. Making this choice 
will simplify calculations considerably, since one can now replace everywhere above k with 
— k and apply symmetries much more efficiently. With this simplification one finds 

ff = o 

g-l = -4A^ [{A^ + B^) - • cos 0] \s\'^ (132) 



g • g = 16^^ sin^ (f) [(C2)2 - (C^ • C 



s\ 



and after some algebra the determinant in Eq. (|129|) takes the embarrassingly simple form 
2^^^|s|^. The logarithm of Eq. ( [T26| ) in the energy expression then becomes 



lny^2M8|s|4 = ln|s|2 + ln4A^, (133) 

where the second term is obviously the extra constant term, whose form depends on the 
choice of the blockage. This has to be retained in order to have convergence in the Mat- 
subara summation, and the final form for the coupling energy is 

i^pinhoie = ]^AohvYNYTTkBTY^ j ^|k-n| [lu | s(k, e„) | ^ + ln(4A4)] . (134) 
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B Important constants 



All calculations presented in the text were done at vapor pressure, which is the situation 
that was present in the experiment of Ref. |13|. Here I just tabulate the values of some 



important constants in the SI units, and give some definitions. First of all we have 

ms = 5.00812 • lO'^^ kg 
h = 1.054573 • 10"^'' Js 
kB = 1.380658 • 10^2^ J/K 

where ms is the mass of a ^He atom, h is Planck's constant divided by 27r and is 
Botzmann's constant. Secondly, for the Fermi wavenumber and Fermi velocity we have the 
relation kp = {m* /h)vF, where m* is the effective mass of a ^He quasiparticle. At vapor 
pressure we have the values 

m* = 2.8m3 
vp = 59.03 m/s 
ksTc = 1.28 • 10"^^ J, 

where Tc is the critical temperature for pure '^He. The density of states at the Fermi surface 
-^(■^k = 0) has the expression Np = = ^^-i^i , which gives, at vapor pressure, the 

following energy and current units: 

Im^VFNpkBTc = 3.7946 kg/m^ s 
HvFNpkBTc = 3.9952 • 10~* J/m^ 
2m3/h = 9.4979 • 10^ kg/Js. 

For easier conversion between Joules and electron volts, I also add the relation 

1 eV = 1.6021773 • 10"^^ J. 
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